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Abstract 


The download on the wing produced by the rotor wake of a tilt rotor vehicle in 
hover is of major concern because of its severe impact on payload-carrying capability. 
In a concerted effort to understand the fundamental fluid dynamics that causes 
this download, and to help find ways to reduce it, computational fluid dynamics is 
employed to study this problem. The thin-layer Navier-Stokes equations are used to 
describe the flow, and an implicit, finite difference numerical algorithm is the method 
of solution. This report is a summary of the methodology developed to analyze the 
tilt rotor flowfield. Included are discussions of computations of an airfoil and wing in 
freestream flows at -90 degrees, a rotor alone, and wing/rotor interaction in two and 
three dimensions. Preliminary results demonstrate the feasibility and great potential 
of the present approach. Recommendations are made for both near-term and far- 
term improvements to the method. 



Contents 


List of Figures iii 

List of Symbols vii 

1 Introduction 1 

1.1 Motivation 1 

1.2 Previous Work 3 

1.2.1 Experimental Work 3 

1.2.2 Theoretical Work 6 

1.3 Current Approach 10 

2 Governing Fluid Dynamic Equations and Solution Method 11 

2.1 General Comments 11 

2.2 Governing Equations 11 

2.3 Turbulence Model 16 

2.4 Numerical Algorithm 17 

2.5 Additional Features 23 

3 Grid Generation 25 

3.1 General Philosophy 25 

3.2 Elliptic Grid Generation 27 

4 Boundary Conditions 30 






4.1 


General Remarks 


30 


4.2 Non-Rotor Boundary Conditions 31 

4.3 Modeling the Rotor 34 

4.3.1 Approach 34 

4.3.2 Combined Momentum Conservation - Blade Element Analysis 35 

4.3.3 Grid Used for Modeling the Rotor 41 

4.4 Wall Jet 43 

4.5 Initial Conditions 44 

5 Discussion of Results 45 

5.1 Preliminary Comments 45 

5.2 Freestream Calculations (No Rotor) 46 

5.2.1 Results in Two Dimensions 46 

5.2.2 Results in Three Dimensions 51 

5.3 Rotor Alone in Two and Three Dimensions 56 

5.4 Wing/Rotor Interaction 65 

5.4.1 Results in Two Dimensions 65 

\ 

5.4.2 Results in Three Dimensions 71 

5.5 Wall Jet in Two Dimensions 78 

6 Conclusions and Recommendations 82 

6.1 Conclusions 82 

6.2 Recommendations S3 

6.2.1 Near Term 83 

6.2.2 Far Term 84 


Bibliography 


86 




List of Figures 

1.1 Sketch of the V-22 in hover, showing the main flow features 4 

1.2 Effect of flap angle on download (Ref. [2]) 5 

1.3 Typical surface pressure distributions on a circulation control airfoil 

in hover mode with and without blowing 7 

3.1 Exponential grid point stretching applied to an arbitrary curve. ... 26 

4.1 Relative velocity and force diagram at a rotor blade element 36 

4.2 Sample output from the momentum theory/blade element analysis for 

a 7-ft diameter, 3-bladed rotor at 2090 RPM and with Ct = 0.0065. . 41 

4.3 One-point overlap grid for wall jet boundary conditions 43 

5.1 Views of a 2-D O-grid around the V-22 airfoil 47 

5.2 Velocity vector plots about the V-22 airfoil at two different points in 

time, in a = 0.2, a = —90° flow 48 

5.3 Close-up view of the velocity vectors around the V-22 airfoil leading 

edge, in a = 0.2, a = —90° flow 50 

5.4 Mach number contours about the V-22 airfoil, in a = 0.2, a = 

-90° flow 50 

5.5 Surface pressure coefficient distribution on the V-22 airfoil, in a = 

0.2, q = -90° flow 51 

5.6 Views of a 3-D grid around the finite-span wing 53 




5.7 Velocity vectors in a vertical plane running spanwise through the wing 

mid-chord, in a M = 0.2, a = —90° flow 54 

5.8 Mach number contours in a vertical plane running spanwise through 

the wing mid-chord, in a = 0.2, a = —90° flow 54 

5.9 Velocity vectors just above the wing and beyond the wing in the x — y 

plane, in a Moo = 0.2, a = —90° flow 55 

5.10 Typical chordwise surface pressure coefficient distributions at various 

spanwise stations along the wing, in a M = 0.2 , a = —90° flow. ... 56 

5.11 Spanwise upper surface pressure coefficient distributions at seven dif- 
ferent chordwise locations along the wing, in a M = 0.2, a = —90° 

flow 57 

5.12 Particle traces released in a flow through a 2-D actuator disk, where 

Ap = O.Olpoo 59 

5.13 Pressure contours for a flow through a 2-D actuator disk, where Ap = 

O.Olpoo 60 

5.14 Velocity vectors in a vertical plane through the 3-D actuator disk, 

where Ap = O.OSpoo 61 

5.15 Particle traces in a vertical plane through the 3-D actuator disk, where 

Ap = 0.05poo 62 

5.16 Mach number contours in a vertical plane through the 3-D actuator 

disk, where Ap = 0.05poc 63 

5.17 Mach number contours in the rotor plane, where Ap = 0.05poo 64 

5.18 Velocity vectors in a horizontal plane immediately downstream of 
the rotor plane showing the effect of swirl, where Ap = 0.05poo and 

M swirl — 0-05 64 

5.19 Particle traces projected on a vertical plane through the actuator disk 

showing the effect of swirl, where Ap = 0.05poo and M 3 wiri = 0.05. . . 65 


IV 




5.20 Farfield view of the 2-D, two-zone grid used for computations of wing/ rotor 

interaction 66 

5.21 Nearfield view of the 2-D, two-zone grid used for computations of 

wing/rotor interaction 67 

5.22 Farfield view of the velocity vectors for 2-D wing/rotor interaction, 

where A p = O.OSpoo and the rotor diameter is 4.78 68 

5.23 Nearfield view of the velocity vectors for 2-D wing/rotor interaction, 

where A p = 0.05poo and the rotor diameter is 4.78 69 

5.24 Mach number contours for 2-D wing/rotor interaction, where A p = 

0.05poo and the rotor diameter is 4.78 70 

5.25 Pressure contours for 2-D wing/rotor interaction, where Ap = O.OSpoo 

and the rotor diameter is 4.78 71 

5.26 Pressure distribution on the surface of the airfoil for 2-D wing/rotor 

interaction, where A p = O.OSpoo and the rotor diameter is 4.78. ... 72 

5.27 Farfield, cutaway view of the 3-D, two- zone grid used for computations 

of wing/rotor interaction 72 

5.28 Nearfield, cutaway view of the 3-D, two-zone grid used for computa- 
tions of wing/rotor interaction 73 

5.29 View in the plane of the rotor of the 3-D, two-zone grid used for 

computations of wing/rotor interaction 74 

5.30 View of a vertical x — z plane of the grid outboard of the wing tip of 

the 3-D, two-zone grid used for computations of wing/rotor interaction. 74 

5.31 Mach number contours in a vertical plane through the mid semi-span 

location for 3-D wing/rotor interaction, where Ap = O.Olpoo and the 
rotor diameter is 4.78 75 


v 




5.32 Pressure contours in a vertical plane through the mid semi-span loca- 

tion for 3-D wing/rotor interaction, where A p = O.Olpoo and the rotor 
diameter is 4.78 76 

5.33 Velocity vectors in a vertical, spanwise plane through the wing mid- 

chord for 3-D wing/rotor interaction, where A p = O.Olpoo and the 
rotor diameter is 4.78 77 

5.34 Velocity vectors around the wing tip in a vertical, spanwise plane 
through the wing mid-chord for 3-D wing/rotor interaction, where 

Ap = O.Olpoo and the rotor diameter is 4.78 78 

5.35 Mach number contours in a vertical, spanwise plane through the wing 
mid-chord for 3-D wing/rotor interaction, where Ap = O.Olpoo and 

the rotor diameter is 4.78 79 

5.36 Mach number contours in a vertical x — z plane beyond the wing tip 

for 3-D wing/rotor interaction, where Ap = O.Olpoo and the rotor 
diameter is 4.78 - - - 79 

5.37 Velocity vectors for a 2-D wing/rotor interaction case with tangential 

wall blowing, where the plenum pressure is l.OSpoo, Ap = 0.05poo 
across the rotor, and = 0.001 80 

5.38 Close-up view of velocity vectors for a 2-D wing/rotor interaction case 
with tangential wall blowing, where the plenum pressure is l.OSpoo, 

Ap = 0.03poo across the rotor, and Moo — 0*2 81 


vi 




List of Symbols 


a 

A 

b 

B 

c 


C P 

C 

Cp 

D$, D r}} Dq 


e* 


E 

F 

G 

G v 

I 

J 

k 

Lz,L v , L { 
M 

Moo 

P 

Pr 

Q 


speed of sound 

rotor disk area, or flux Jacobian matrix in the ^-direction 
wing span 

flux Jacobian matrix in the //-direction 
wing chord 

specific heat at constant pressure 
flux Jacobian matrix in the (-direction 
pressure coefficient 

diagonal matrices associated with the LU-ADI algorithm 

total energy per unit volume 

internal energy per unit mass 

inviscid flux vector in the (-direction 

inviscid flux vector in the //-direction 

inviscid flux vector in the (-direction' 

viscous flux vector in the (-direction 

identity matrix 

transformation Jacobian 

coefficient of thermal conductivity 

lower bidiagonal matrices associated with the LU-ADI algorithm 

viscous flux Jacobian matrix in the (-direction 

freestream Mach number 

pressure 

Prandtl number 

vector of conserved quantities 


vn 




R 

Re 

t 

T(,T„T ( 

W, V, w 

u t ,u, f u< 

u,v,w 

x,y,z 


radial location 
rotor radius 
Reynolds number 
time 

similarity transformation matrices 
Cartesian velocity components 

upper bidiagonal matrices associated with the LU-ADI algorithm 
contravariant velocity components 
Cartesian coordinates 


Or 


7 

a 4 ,a„,a c 
W, v c 

e 

A< 

A, 

Ac 


P 


u> 


angle of incidence 
ratio of specific heats 
central difference operators 
forward difference operators 
backward difference operators 
artificial dissipation coefficient 

chordwise, spamvise, and normal coordinates in body conforming system 
diagonal matrix of eigenvalues associated with the flux Jacobian matrix A 

diagonal matrix of eigenvalues associated with the flux Jacobian matrix B 

diagonal matrix of eigenvalues associated with the flux Jacobian matrix C 

coefficient of viscosity 
density 

angular frequency of rotation 


viii 




Chapter 1 
Introduction 


1.1 Motivation 

The tilt rotor aircraft is a flight vehicle which combines the vertical takeoff and 
landing capability of the helicopter with the efficient high-speed cruise performance 
of conventional fixed-wing aircraft. The tilt rotor vehicle has a rotor located at each 
wing tip. The rotors can be rotated to permit conversion from the helicopter mode 
to the airplane mode and vice versa, the wing remaining fixed. 

The concept was first proposed by Bell Helicopter engineers during World War II, 
and it evolved into a first prototype in 1955, designated the XV-3 [1]. In 1977, the 
NASA /Army /Bell XV- 15, a 13,000 lb experimental tilt rotor aircraft, flew for the 
first time in a research program that continues to today. The usefulness of the tilt 
rotor aircraft is evidenced in the recent development of the V-22 Osprey for the 
U.S. Armed Forces by a Bell Helicopter Textron/Boeing Vertol team. The V-22 is 
a multi-service, multi-mission tilt rotor aircraft. It has a vertical take-off weight of 
47,500 lb and is capable of transporting up to 40 passengers. 

The unique features of the tilt rotor can also be exploited as a civil transport in 
the city-center to city-center commuter market or as a feeder to hub airports. The 
need for such a mode of transport will certainly increase as community real estate 
prices continue to increase, making new airport construction prohibitively expensive, 
driving new airport locations further away from large population densities. 
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The tilt rotor concept offers several considerable advantages over the rival tilt 
wing concept (in which the rotors and wing both rotate in the transition from heli- 
copter to airplane mode and back). Wing tilt requires additional mechanical com- 
plexity and resulting increased structural weight to support the higher concentrated 
wing/fuselage junction loads. The tilt rotor concept also offers much greater con- 
trollability while hovering in gusty wind or cross-wind conditions. 

A major limitation of the tilt rotor configuration, however, is the aerodynamic 
download imposed on the wing by the rotor flowfield when hovering. Because the 
wing is fixed, the rotor flow, in hover, hits the wing at 90 degrees. The download 
force on the wing has been measured and can be as large as 10 - 15 percent of 
the total rotor thrust [2,3]. Assuming the payload-carrying capability to be about 
25% of gross take-off weight, complete elimination of the download could increase 
the effective payload by over 50%. The need for a thorough understanding of, and 
the eventual reduction of, wing download, then, is the major impetus driving this 
theoretical study on tilt rotor flowfields. 

Study of the three-dimensional wing/rotor interaction, in hover, therefore, is the 
primary focus of this work. The flowfield about a tilt rotor configuration is very 
complicated. The rotor, typically located less than a wing chord above the wing, 
induces a flow which is closely coupled to the flow about the wing. There exists a 
large region of nearly-stagnated flow on the wing upper surface and a large region of 
unsteady, turbulent, separated flow below the wing. The flow over the wing upper 
surface is highly three-dimensional. Near the wing tip, the flow over the wing is 
essentially chordwise. Further inboard, the flow becomes increasingly spanwise. Due 
to symmetry of a hovering tilt rotor, the spanwise flow from both wings meets at the 
vehicle centerline and is redirected upwards. Some of this rising column of air is re- 
ingested by the rotor thus creating a large-scale recirculation pattern which reduces 
rotor performance. This flow pattern has been termed the “fountain effect” . Refer 



to Figure 1.1 for a simplified sketch of the main flow features about a V-22 in hover. 


1.2 Previous Work 

1.2.1 Experimental Work 

Flight test of the XV- 15 [1,4] has yielded quantitative estimates of hover performance 
including download as a function of flap angle. Figure 1.2 shows the download (DL) 
normalized by the rotor thrust (T) plotted as a function of flap angle. Flight test 
results are compared to estimates from a semi-empirical theory by Felker and Light 
[2] . To better study the nature of the tilt rotor flowfield itself, the flexibility and 
control offered by wind tunnel testing is required. McCroskey et al. [5] measured the 
drag of two-dimensional models of the XV- 15 airfoil with various flap and leading 
edge configurations. They found that the drag on the airfoil in a freestream flow at 
-90 degrees was very sensitive to not only flap angle but also the surface curvature 
distribution on the upper surface near the leading edge. Lowering flap reduces the 
frontal area thereb\ r reducing the download. A flat plate has a 2-D drag coefficient 
about twice that of a circular cylinder. Increasing airfoil thickness and camber, then, 
which tend to make the airfoil less like a flat plate and more like a circular cylinder 
or ellipse, reduces the download. 

Maisel et al. [6] continued this 2-D experimental effort by examining the effects 
of several different airfoil, flap, and leading edge configurations on the download. 

Boeing has tested a powered tilt rotor model whose basic geometry was that of 
a 0.15 scale V-22 Osprey. Some results of this test are reported in reference [7], 

Results from model tests of tilt rotor hover and wing/rotor interaction by NASA 
have been reported in references [1,2,3]. Large-scale tests involving a 0.658 scale 
V-22 rotor and wing as well as small-scale tests of a 0.16 scale S-76 helicopter rotor 
and a model wing were carried out at the Outdoor Aerodynamic Research Facility 
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Figure 1.1: Sketch of the V-22 in hover, showing the main flow features. 
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Figure 1.2: Effect of flap angle on download (Ref. [2]). 

(OARF) at NASA Ames Research Center. These tests provided flow visualization 
of the complicated flowfield as well as wing pressure measurements and estimates 
of wing download and rotor performance. A clearer picture of the flowfield was 
obtained. The gradual transition from near-2-D flow on the outboard section of 
the wing to near-spanwise flow on the inboard region, and the fountain effect were 
observed. 

Because the static pressure in separated flow is generally somewhat less than 
ambient pressure, a download suction force on the wing lower surface contributes to 
the total download. By energizing the boundary layer using tangential blowing on 
the upper surface near the leading edge (where the flow is close to separating) the 
separation location can be further extended around the leading edge to the lower 
surface thereby reducing the chordwise extent of the of the separated flow region 
below the wing and significantly reducing the download. This effect was observed 
in the experimental work described in [2] using the 0.16 scale model and wing with 
elliptic airfoil section with blowing slots. Figure 1.3 shows typical surface pressure 
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distributions on a circulation control airfoil with and without blowing in a rotor- 
induced -90 degree slipstream.. The slots were located at 3% and 97% chord, the 
former blowing towards the leading edge, the latter towards the trailing edge. It 
was found that particularly the leading edge slot was very effective in controlling the 
location of separation and in reducing the download. 

A more detailed set of tests similar to those described above, scheduled for the 
Fall of 1989, will soon provide even more information. A 0.18 scale V-22 rotor 
and both a V-22 model wing and an XV-15 wing with blowing slots will be tested. 
This experimental research effort will provide a means of direct comparison with the 
computational results obtained in this theoretical study, whose preliminary results 
are reported here. 

1.2.2 Theoretical Work 

The theoretical prediction of tilt rotor flows and of download is very difficult because 
of the complexity of the flowfield. The rotor wake is an unsteady, three-dimensional 
flow with regions of concentrated vorticity (rotor tip vortices). The rotor blades 
themselves see highly compressible flows, and beneath the wing there exists a large 
region of unsteady, turbulent, separated flow. Accurate, simultaneous prediction of 
all these flow features still lies beyond the state of the art. Simplifications must be 
made to render the problem tractable. 

References [8,9] describe the application of a low-order panel method to the tilt 
rotor problem. The rotor was modeled using source singularities and the rotor wake 
was represented by a time-averaged cylindrical vortex sheath comprised of quadrilat- 
eral doublet panels. A blade element model of the rotor feeds time-averaged loading 
as a function of radial and azimuthal location to the panel code which contains a 
model of the wing. The wing was modeled simply as a cambered plate using a lattice 
of doublet singularities. Many of the overall flow features were predicted using this 
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Figure 1.3: Typical surface pressure distributions on a circulation control airfoil in 
hover mode with and without blowing. 
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method. Quantitative results, however, because of the nature of the equations solved 
(Laplace’s equation) must be viewed with caution as separated flows cannot be accu- 
rately predicted with this formulation. As found in the experimental work described 
in reference [5], the separation location is very sensitive to leading edge curvature 
and thickness. These important effects cannot be accurately predicted using a panel 
method. Download, which is dependent on viscous effects, can only be accurately 
predicted using an analysis which incorporates the effect of viscosity. 

References [5,10] describe discrete- vortex seeding methods to calculate the un- 
steady, 2-D flow around an airfoil at an angle of incidence of -90 degrees. In reference 
[5], the wing is immersed in a freestream flow. In reference [10], to study rotor/wing 
interaction, a rotor is modeled using constant strength doublet panels which induce 
a normal velocity distribution. Since no integral boundary layer calculation was 
coupled with the potential flow calculation, boundary layer growth and separation 
location were not predicted. Separation location was specified and a uniform base 
pressure on the wing lower surface was assumed. The methods predicted the upper 
surface pressures fairly well but obviously were incapable of accurately calculating 
the lower surface pressure. Since the flow over the tilt rotor wing is highly three- 
dimensional, two-dimensional analyses such as these are of limited usefulness. 

Recent developments in numerical algorithms and in computer speed and memory 
capability have made tractable the solution of the Euler and Navier- Stokes equations 
for increasingly complicated flows. However, only very few attempts at the calcu- 
lation of rotor flows have so far been made using the latest computational fluid 
dynamics (CFD) techniques. 

Quasi-steady solutions of a 2-bladed rotor in hover have been obtained by Srini- 
vasan and McCroskey [11] using a flux-split, approximately-factored, implicit algo- 
rithm to solve the unsteady, thin-layer Navier-Stokes equations. The computation 
of the concentrated vortices shed from the rotor tips is affected adversely by nu- 
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merical diffusion due to insufficient grid density beneath the rotor. Therefore, the 
effect of the rotor wake, and in particular, the induced downwash, must be estimated 
empirically in order to obtain reasonable results for lifting cases. 

The same authors in reference [12] presented their results of airfoil/vortex inter- 
action in two-dimensions. They describe a method where the structure of the vortex 
is prescribed, but its path in space is allowed to develop as part of the solution. 
Using this approach, they were able to simulate the effect of a shed tip vortex on a 
section of the following rotor blade. 

Roberts [13] combined an unsteady Euler solver with a wake model for a two- 
bladed hovering rotor. The bound circulation distribution along the span of each 
rotor blade was determined from the Euler solution and used to set the strength of 
the wake vortices. 

McCroskey and Baeder [14] estimated that in order to calculate two revolutions of 
a two-bladed rotor above a simple fuselage using a typical, implicit, thin-layer Navier- 
Stokes code with algebraic eddy- viscosity modeling of turbulence, a 100 megaflop 
computer would require 40 CPU hours and 30 million words of memory ( or 4 hours 
of CPU for a one gigaflop machine). Routine calculations of 3-D rotorcraft flows 
including detailed modeling of the rotor blades will remain elusive for some time. 

Rajagopalan and Mathur [15] modeled a three-dimensional rotor in forward flight 
using a distribution of momentum sources added to the steady, incompressible, lami- 
nar Navier-Stokes equations. Rotor geometry and blade sectional aerodynamic char- 
acteristics were incorporated into the evaluation of the source terms. Their results 
represent a time-averaged solution. In complexity, this method lies between an ac- 
tuator disk representation and a CFD computation of the individual blades. 
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1.3 Current Approach 


Despite the research efforts of the past several years, gaps in our fundamental un- 
derstanding of the tilt rotor fiowfield remain. It is the objective of this current work 
to model the tilt rotor fiowfield as accurately as present computational tools allow. 
Computer resources and solution algorithms have evolved, during the past several 
years, to the point where the solution of the thin-layer Navier-Stokes equations about 
a simplified tilt rotor configuration comprised of wing and rotor model is technically 
feasible. Since the modeling of the flow about the individual rotor blades remains 
such a complex task and tremendous computational drain, and since the primary in- 
terest here is wing download prediction and not detailed rotor simulation, the rotor 
in this study, is modeled as an actuator disk. The loads from the blades are aver- 
aged over the entire rotor disk area. The viscous flow around the airfoil, however, is 
accurately predicted by defining a suitable grid. A multiple zone approach is used 
to incorporate the effects of the rotor. Empirical data or, in this case, a momen- 
tum theory/blade element analysis is used to estimate the average radial distribution 
of axial and swirl velocities, and the pressure rise across the rotor disk (which are 
required for the CFD calculation). 

The theory and numerical method, grid generation, development of boundary 
conditions, and preliminary results obtained, will all be discussed in greater detail 
in the chapters to follow. 
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Chapter 2 


Governing Fluid Dynamic 
Equations and Solution Method 

2.1 General Comments 

The work described in this report involves the numerical solution of the three- 
dimensional Navier-Stokes equations. It is assumed that there are no body forces 
(eg. gravity effects are not important) and there is no heat addition or removal. The 
thin-layer approximation, described later, is applied. It assumes that the viscous 
forces are confined to a small region near the wing surface. The solution method, in 
the form in which it is applied here, was formulated by Fujii and Obayashi [16,17,18]. 
It is suitable for subsonic, transonic, and supersonic flow calculations. Although the 
solution method is implicit, the boundary conditions are applied explicitly. This 
allows easier application of the code to a wide variety of geometries. 

The general formulation is outlined in the following sections. Greater detail can 
be found in references [16,17,18,19]. 

2.2 Governing Equations 

The Navier-Stokes equations are the most basic continuum-representation of fluid 
dynamic flows. To ensure proper shock capturing (in transonic flow, for example), 
i.e. to accurately predict shock strength and location, the three-dimensional Navier- 
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Stokes equations are written and solved in strong conservation law form where the 
dependent variables are expressed in the form of spatial gradients. These equations 
can be seen in several different references. See, for example, reference [19]. To 
convert them to a more useful form for computational purposes, the equations are 
manipulated somewhat. This is described below. 

The equations are non-dimensionalized by dividing the density p by the freestream 
density p^, the velocity components u, v , and w by the freestream speed of sound 
aoo, and the total energy per unit volume e by p<x>alo- Conforming to the normal 
convention, u is in the downstream (x) direction, v is in the spanwuse (j/) direction, 
and w is in the vertical, upwards (c) direction. The coefficient of viscosity is normal- 
ized by and the lime is normalized by c/a ^ w-here c is the wing chord. Applying 
this non-dimensionalization to the Navier-Stokes equations results in a term con- 
taining the expression {pooa<x>c)l P<x>- introducing the definition of Mach number, 
Moo = Uoo/aoo where u ^ = y/v? + u 2 + v? , then it is seen that this expression 
equals Re/Moo where Re = (poo^ooc)//ioo- 

Next, the Navier-Stokes equations are transformed from Cartesian coordinates 
to general curvilinear coordinates. This makes tl\e formulation independent of the 
body geometry thereby easing the specification of the boundary conditions, ft also 
allows for straight-forward application of the thin-layer assumption later. In addition, 
since the physical domain is transformed into the computational domain which is a 
rectangular box with a uniformly-spaced mesh, then standard differencing schemes 
can be used for the derivatives. The coordinate transformation is defined by: 

T = t 

Z = fiixiy,*,*) 

V = I i(x,y,z,t) (2.1) 

C = C (x,y,z,t) 
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where t represents the time independent variable. The airfoil surface is transformed 
to the ^-direction, the spanwise direction is transformed to 77 , and ( is roughly normal 
to the wing surface. The governing equations must, correspondingly, be transformed 
from the physical domain (x,y,z) to the computational domain (£, Details 

of this transformation procedure can be found in references [20,21]. By writing 
the transformation in terms of spatial derivatives and applying the chain rule, a 
transformation Jacobian, J, and several identities called metrics can be defined as 
follows: 


J — 1 / det 


X v X C 
% Vn Vi 


( 2 . 2 ) 


where x ^ — dx/d etc. i.e. all the above matrix elements area partial derivatives. 
They are evaluated numerically using second-order accurate differences. The metrics 


are: 


— J ivv z c 

Zy = J ( x ( z v ~ x n z <)i 

iz = J c ynh 
r)x = J (y<Zz -VtZ (), 


Cx = J -yv z z) 


C y = J ( x v z i - x t z v) 


Cz — J { x iyr) x 7)yt) 

£,t = x t£x 2/r£y Z r£z 


Vy — J ( X £ ~C X <i Z 0'' Vr f]y ZtVz 

Vz = J { X CVt Ct = ~ x tCx ~ VrCy z tCz 


(2.3) 


Note that for stationary grids (no body motion), the metric time derivative terms 
are zero. 

Generally, the thin-layer approximation is applied to the Navier-Stokes equations 
to reduce the computational effort (particularly in three dimensions) to a manageable 
level. For the high Reynolds number flows which are typical of practical aerodynamic 
problems, the viscous effects are confined to a small region near the body surface 
and in the wake. Computer memory limitations usually necessitate concentrating 
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the available grid points near the surface of interest. This results in fine grid spac- 
ing normal and near to the surface and coarse grid spacing tangential to the surface. 
With this type of grid, even if the full Navier-Stokes equations were programmed, the 
viscous terms associated with the flux vectors along the body would not be resolved, 
and for most cases of interest, these terms are negligible anyway. Therefore, it is jus- 
tifiable to eliminate from the calculation the viscous fluxes in the directions parallel 
to the surface, i.e. the £ (chordwise) direction and the tj (spanwise) direction. This is 
easy to do with the equations already transformed to the body-fitted computational 
domain. The thin-layer approximation is similar in philosophy but somewhat differ- 
ent than the assumptions in boundary layer theory. Less restrictive than boundary 
layer theory, the thin-layer approximation retains the normal momentum equation 
and allows pressure variation through the boundary layer. 

The thin-layer approximation is less valid for low Reynolds number flows and 
in regions of large flow separation. This then is one of the major limitations of 
this method for studying the 3-D tilt rotor flowfield. Even if the full Navier-Stokes 
were solved, the limitations of the turbulence model in regions of extensive separation 
would contribute to inaccuracies in the computed flowfield in this region. The current 
approach, however, makes the calculation tractable and far superior to any method 
hitherto applied to this problem. The full 3-D Navier-Stokes solution is beyond the 
means of current practically-available computer power. 

Applying the thin-layer approximation, then, the non-dimensional, 3-D Navier- 
Stokes equations in transformed curvilinear coordinates become: 

dQ dE 8F dG M ao dG v 

dr + di + di] + d( ~ Re d( ( ' } 

The Q vector contains the transformed conservative flow variables 

P 

pu 

Q = J~' pv 

pw 
e 
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Note that the elements of the Q vector, as well as all flow variables referred to in 
subsequent discussions in this chapter, are non-dimensional quantities. 

E , F, and G are vectors that contain the inviscid “Euler” terms: 


E= J~ l 


F — J 


-l 


P U 

puU + i x p 
PVU + ZyP 
pwU + £ z p 
L U (e + p) - ( t p J 

pV 

pilV + Tj x p 
pvV + T)yP 

pwV -t- rj z p 
V (e + p) -7 pp 


( 2 . 6 ) 


(2.7) 


p\V 

puW + ( x p 

G = J -1 | pvW + CyP | (2.8) 

pwW + CzP 
W (e + p) - CtP 

where p is the static pressure and £/, V, and IT are the contravariant velocity com- 
ponents that appear as a result of the coordinate transformation. They are defined 
below: 


U = £t + £xU + iy v + CzW 

V = Tjt + Tj x U + T}yV + r) z w (2.9) 

V = (t + + Cy 17 + (z w 


The vector G v contains the viscous terms retained after the application of the thin- 


layer approximation. It is given by 


where 


Gy - J 


Mo 


Re 


0 

Cx^rx T Cy T xy T CzTxz 

Cx T yx T Cy T yy T Cz^yz 

CxTzx T Cy^zy T Cz^zz 

CxPx + (yPy + (zflz 


2 

r rT — "p(u r + v y + w z ) + 2fxu T 

o 


( 2 . 10 ) 
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T yy = ~~H{U X + Vy + W Z ) + 2fiVy 

2 

Tzz = + Vy + Wz) + ^HW Z 

T X y Ty X /i (liy T Vj; ) 

Txz = T zx = H ( U 2 + W x ) 

Tyz = T zy = H (V z + W y ) 


fix 

fiy 

fiz 


jj, 

~~~ O x €-i “I - 11T XX T VT X y T IVT X 

Pr 


2t_ 

Pr 


dy C{ “j - UTy X “I" VTyy T IV Tj 


yz 


1 1 

—d z e t + UT ZX + VT Z y + WT Z 

Pr 


The internal energy per unit mass, is 


e ( u 2 + v 2 + w 2 ) 

6 ' = ~p 2 


(2.11) 


The Prandtl number Pr is defined as Pr — c p p/k where c p is the specific heat at 
constant pressure and k is the coefficient of thermal conductivity. Also, 7 is the ratio 
of specific heats. For air, at temperatures and pressures of interest here, Pr = 0.72 
and 7 = 1.4. Pressure is related to the conservative flow variables through the 
equation of state for a perfect gas 


P = (7 - 1 


e 



( 2 . 12 ) 


To evaluate the spatial derivatives of the Cartesian velocity components in Eq. 2.11, 
the chain rule is applied. For example, 


u x = £ X U(. + T) x Ur, + Cr«< 


2.3 Turbulence Model 


In order to resolve the turbulent eddies which contribute to a flow’s turbulence level, 
an extremely tiny grid spacing and a huge number of points would be required. 
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For the solution of practical problems of interest to be feasible, a turbulence model 
is generally introduced which greatly reduces the number of grid points required 
in viscous regions. Based on a combination of theory and empiricism, turbulence 
models, although far from precise predictors of turbulence, do provide a means of 
modeling a real flow with a turbulent boundary layer. 

For the purposes of this preliminary work, no turbulence model was incorporated. 
Only laminar flow calculations were performed because the main focus was on mod- 
ification of the Navier-Stokes code for the tilt rotor problem, and on generation of 
suitable grids and boundary conditions. 

A turbulence model proposed by Baldwin and Lomax [22] for conventional bound- 
ary layers and another for turbulent wall jets on curved surfaces proposed by Roberts 
[23] will be incorporated in the near future. 


2.4 Numerical Algorithm 


The solution technique employs an implicit, approximately factored, non-iterative 
method first developed by Beam and Warming [24]. Explicit methods suffer the 
disadvantage of having a severe restriction on time step size in order to maintain 
stability. This is particularly acute for Navier-Stokes solutions where, because of 
the relatively small scales associated with resolving the boundary layer, the partial 
differential equations are very stiff. Often, the steady-state solution is of principal 
interest, so being able to use large time steps to accelerate the rate of convergence is 
very important. Implicit methods are stable for relatively large time steps even for 
highly nonlinear equations such as the Navier-Stokes equations. 

Applying the first-order accurate (in time) implicit Euler scheme to Eq. 2.4 results 


in 


Q n+ 1 - Q n + to 


0E n+ 1 dF n+l <9G n+1 Mqo dG” +1 


+ 


+ 


= 0 


(2.13) 


' dr, ' d( Re d( 
where n + 1 is the time at which Q is desired, n is the previous time level at which 
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Q is known everywhere, and Q n = Q(nAi ) . 

Since the flux vectors E, F, G , and G v are nonlinear functions of Q, then Eq. 2.13 
is nonlinear in Q n+1 . In order for the method to be non-iterative, the nonlinear terms 
are linearized in time about Q n by a Taylor series expansion such that 

E n+ 1 = E n + A n AQ n + O (At 2 ) 

F n+i _ + B n AQ n + O (At 2 ) (2.14) 

G n +i = c n + C n AQ n + O (At 2 ) 

G'” +1 = G n v + M n AQ n + O (At 2 ) 

Note that AQ n = Q n+1 — Q n . Also, A n , B n , C”, and M" are the flux Jacobians. 
Expressions for these matrices can be found in reference [25]. 

In the Beam and Warming method, the alternating direction implicit (ADI) al- 
gorithm replaces the inversion of one huge matrix — which would be prohibitively 
expensive to compute — to the inversions of three block tridiagonal matrices, one 
for each direction. Efficient block tridiagonal inversion routines exist, making this 
algorithm a viable solution technique. 

A linear constant coefficient Fourier stability analysis (assumes periodic boundary 
conditions) for the three-dimensional model wave equation shows a mild, uncondi- 
tional instability for the Beam- Warming factored algorithm [26]. A small amount 
of artificial dissipation (also called “smoothing”) is required to render the scheme 
stable. The most common procedure is to add fourth-order artificial dissipation to 
the explicit right hand side of the equation (see Eq. 2.15 below) and second-order 
dissipation to each of the three implicit one-dimensional block operators on the left 
hand side. Second-order implicit dissipation is used so as to keep the block implicit 
operators tridiagonal. 

Applying the Beam- Warming scheme to Eq. 2.13 using the linearized expressions 
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for the flux vectors, Eqs. 2.14, the solution algorithm can be written as 

[/ + AtS^A 71 - [l + M6 t ,B a -Att I J- l V v A r ,j] 

M 

I + At S ( C n - At - At e I J~ 1 V r) A r ,J A Q n 

He 

r M 1 

= -At «<£" + 6,F' + 6<G“ - (215) 

-At e E J-' [C V< At ) ! + ( V„A„) J + (V< A<) J ] J Q " 

In the above notation, 

[/ + At <^A n ] A Q n = A Q n + A t6 t {A n AQ n ) 

Note that / is the identity matrix, <5 is a central difference operator, and V and A 
are backward and forward difference operators, respectively, tj and ce are constants 
that control the magnitude of artificial dissipation introduced in the scheme on the 
implicit side (left hand side) and on the explicit side (right hand side) of Eq. 2.15, 
respectively. 

The algorithm is said to be in “delta form”. For steady-state solutions ,A Q n — ► 0, 
and the solution is independent of the choice of operators on the left hand side of 
Eq. 2.15. The algorithm is first-order accurate in time. Since second-order spatial 
differences are used on the implicit side, the method is second-order in space. For 
converged, steady-state solutions, when the implicit side (left hand side) of the equa- 
tion approaches zero, and if fourth-order differencing is used on the explicit (right 
hand side terms), then the method becomes fourth-order accurate in space. 

Most of the computational effort involved in an implicit algorithm such as the 
one outlined above is associated with inverting the block tridiagonal matrices in 
each direction. Pulliam and Chaussee [27] suggested a way to diagonalize the blocks 
within the block tridiagonal matrices, thereby transforming them to scalar tridiagonal 
matrices which are much more efficient computationally to invert. Their approach 
is based on the fact that the flux Jacobians A, fi, and C each have real eigenvalues 
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and a complete set of eigenvectors. This means that the Jacobian matrices can be 
diagonalized by similarity transformations as indicated below: 


A £ = r { -MT { , A v = T~ l BT v , A ( = r ( -‘cr ( 

where, for clarity, the superscript n has been dropped from the flux Jacobians. A^, 
A„, and A^ are diagonal matrices containing the eigenvalues of matrices A, B , and 
C, respectively. T%, T v , and are similarity transformation matrices. Expressions 
for the above matrices can be found in references [26,27]. 

For simplicity, looking at only the ^-direction, the Beam- Warming’s ADI implicit 
operator (see Eq. 2.15) can be written in the diagonal form as 


I + AtS^A + A<e/J _1 V € A € J 


= T(_T[ l + A t6t + AtJ-'aVzAtJ 



I + At 6 € A { + At Tf 1 


(2.16) 


Moving and outside of the difference operator 6^ introduces an error which 
renders the method at most first-order accurate in time [27]. For steady-state cal- 
culations, however, where the left hand side of Eq. 2.15 goes to zero as AQ n — ► 0, 
the solution obtained using the diagonal algorithm will be identical to that obtained 
from the original Beam- Warming ADI scheme since the right hand side is identical 
for both methods. 

Fujii, Obayashi, and Kuwahara [16,17,18] introduced a further modification to the 
left hand side operators that reduces the tridiagonal matrix — obtained after the 
diagonalization described above — to a product of a lower and an upper bidiagonal 
matrix, thereby reducing the computational effort even more. This is possible by 
employing the flux vector splitting technique introduced by Steger and Warming in 
reference [28], and by using a diagonally dominant factorization first proposed in 
reference [29]. These modifications are outlined below’. 
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The central differencing of Eq. 2.16 can be decomposed into two one-sided dif- 
ferences by using the flux vector splitting technique. Continuing to look at the 
^-direction operator only, Eq. 2.16 can be re-written as 


I T At 8^A -j- At £/</ * ^7 ^A^J 

* T< [i + V* A+ + A { Af] T- 1 (2.17) 

with 

A? = y (At ± |A C |) ± (2.18) 

where contains all the positive eigenvalues and contains all the negative 
eigenvalues. Note that the positive-moving characteristics (eigenvalues) are backward 
differenced, and the negative- moving characteristics are forward differenced. This is 
required for numerical stability. the inverse of the Jacobian, is evaluated at 

the central point. Using first-order, three-point, upwind differencing (three-point 
differences offer greater numerical stability than two-point differences) for the V^* 
and A^ operators, Eq, 2.17 can be written as 


~ TJr + D £ + U$\ 


where 


h 

Di 

U< 


1 + e ( A i _ A j ) 

8 a _ l._ 

6 Aj+1 6 Aj+2 


where the subscript j is the grid point index in the ^-direction. 
Applying diagonally dominant factorization (refer to [29]), 


L £ + -f — (L^ + Dt) D ^ 1 (D{ + U^) + 0 


(2.19) 


( 2 . 20 ) 


( 2 . 21 ) 
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This can be shown by examining Eqs. 2.18 and 2.20 and noting that D $ is of order 
1 and L{, are of order At. Substituting Eq. 2.21 into Eq. 2.19, Eq. 2.17 becomes 

I -4- At 6^ A -|- At tj 

W + Dz) D^ 1 {Dt + Uz) Tf l (2.22) 

The block tridiagonal system of the left hand side of Eq. 2.16 has been trans- 
formed to matrix multiplications and the product of the lower scalar bidiagonal 
matrix (L^ + D^) and the upper scalar bidiagonal matrix D^ 1 ( D ^ + U^). The same 
analysis is carried out for the ^-direction and a slightly modified analysis is performed 
in the (-direction. 

The viscous flux Jacobian M n is not diagonalizable by the similarity transfor- 
mations and T£~ l and, therefore, the ('-direction implicit operator of Eq. 2.15 
cannot be reduced to diagonal form. To retain the diagonalization in all three direc- 
tions (and not incur the computation penalties associated with not simplifying the 
block tridiagonal operator in the (-direction), the viscous flux Jacobian M n is simply 
neglected. Pulliam [26] , using the unsplit diagonalized method, compared results 
obtained with and without the implicit viscous flux Jacobian. He found that neglect- 
ing this term actually yielded the most efficient steady-state solutions. Pulliam and 
Steger [32] employ this method for steady viscous flows and convection dominated 
unsteady flows. Guruswamy [30,31] , using this same numerical algorithm (ignoring 
the implicit viscous flux Jacobian) for both the Euler and Navier-Stokes equations, 
obtained good results for unsteady aerodynamic and aeroelastic calculations. 

In order to maintain stability of the Fujii /Obayashi method with the thin-layer 
viscous terms retained on the explicit, right hand side of the equation, the diagonal 
matrix of eigenvalues, of the inviscid flux Jacobian matrix ,C n , is modified (see 
references [16,18]) as indicated below (compare with Eq. 2.18, for example): 



where 


tyy/G + Q + Cl 

ik ^ a /• 


where, in the computational domain, A£ = A 7 / = A£ = 1. 
Finally, the present scheme can be summarized as follows, 


+ A) V (A + + a,) k _1 (a, + u*) T ^ T C 

(L< + D<) D? {D c + U c ) A Q n 


= -At 


8^E n + 6 v F n + <5 C G" 


M a 


Re 


a g: 


(2.23) 


—At e E r' [(■ V< Arf + (V.A,) 2 + (VcA c ) 2 ] J Q n 


Analytical expressions for Tr l T v and T ~ 1 and their inverses can be used to reduce 
the computational effort. They are presented in reference [27], 

It is evident from the above equation that the inversion process has been reduced 
to one forward scalar sweep and one backward scalar sweep in each direction, and 
simple matrix multiplications. 


2.5 Additional Features 

To further enhance convergence speeds for steady-state calculations, a space varying 
time step size At can be specified. This modification can be very effective for typical 
grids that have a wide variation of grid spacing. By scaling At with grid spacing, 
a more uniform local Courant number (ratio of local time step to grid cell width 
multiplied by the characteristic velocity) can be maintained throughout the flowfield. 
Since the local transformation Jacobian, J, scales with the inverse of grid cell volume, 
the following has been found to work well ( refer to [32]: 


At = 


A* I,; 

1 + \/j 


where At\ Te j is a fixed, inputted time step. Due to numerical approximations made 
to calculate the metrics, a computed freestream flow may be somewhat in error. To 
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remedy this problem, freestream subtraction is employed [33]. The freestream values 
of all the fluxes are calculated and subtracted from all the computed fluxes. This al- 
lows for recovery of the freestream condition without affecting the basic formulation. 
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Chapter 3 


Grid Generation 

3.1 General Philosophy 

Grid generation is a very important aspect of computational fluid dynamics. The 
grid is the assembly of points at which the numerical solution to the relevant partial 
differential equations is found. To maintain solution accuracy, the grid should possess 
a smooth distribution of points (refer to [34]). Also, the grid (or mesh) and the 
distribution of points must be compatible with the fluid dynamics equations being 
solved and the particular flowfield expected. 

As discussed in Section 2.2, the thin-layer approximation is made so as to limit 
the computer storage requirements to a manageable level. The viscous fluxes normal 
to the surface are dominant. In order to resolve this important contribution, the grid 
spacing must be very fine in a direction normal to the surface. Grid spacing can be 
much coarser along the surface where the far less significant tangential viscous fluxes 
need not be resolved. In addition, grid points must be clustered in regions where 
relatively large flow gradients are anticipated. For example, the wing leading edge 
and trailing edge as well as the mid-chord regions for a = —90° flow require grid 
point concentration. Figure 5.1 is an example of the grid generated for a freestream 
flow computation about the V-22 airfoil at an angle of incidence of -90 degrees. 
There are 47 points in the circumferential direction around the airfoil surface, three 
of which overlap to ease implementation of the periodic boundary conditions (refer 
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Figure 3.1: Exponential grid point stretching applied to an arbitrary curve. 


to Section 4.2 for a further discussion of periodic boundary conditions). There are 
33 points in the direction normal to the airfoil surface. 

To obtain the desired clustering of grid points on a body surface and in the mesh 
interior for algebraically-defined grids, extensive use is made of an exponential type 
of stretching. If a distribution of N points is desired along a curve of specified length 
S, and the arc length between the first two points is specified to be AS, then an 
expression for the total length can be written as (refer to Figure 3.1 ) : 

S = AS + aAS + a 2 AS + a 3 AS + • • • 

N —2 

= A S£a* 

k= 0 


Defining a function / where 

N - 2 

f = s- AS X>* 

k = 0 


then an iterative root finding procedure such as the Newton-Raphson method is 
used to determine the value of a that satisfies / = 0 within a desired tolerance. This 
method can easily be extended to a distribution of points with exponential stretching 
in both directions, i.e. a different AS is specified at both ends of the interval. 

To generate a mesh for a 2-D calculation, five identical 2-D grids were stacked 
parallel to each other. This is required because the Navier-Stokes code employs 
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a fourth-order accurate 5-point spatial numerical differencing. The computational 
results are then referred to as “ pseudo 2-D ” as they are essentially 2-D in char- 
acter although they are generated by a 3-D code. This is made possible by proper 
application of symmetry boundary conditions at both ends of the mesh (refer to 
Section 4.2). 

For actual 3-D cases, 2-D grids are generated at many spanwise stations along 
the wing and beyond the wing tip, and stacked in parallel. 

3.2 Elliptic Grid Generation 

To ensure smooth grid point distributions on the interior of a 2-D mesh, an elliptic 
grid solver can be employed. It is well-suited to the types of grids employed in this 
study. Because of the symmetry of the 2-D flowfield at a given station of the tilt 
rotor wing in hover, an O-grid is more suitable than either a C-grid or an H-grid. 
The O-grid is so named because of the way the grid lines encircle the airfoil. 

The elliptic grid generation scheme was first proposed by Thompson, Thames, 
et al. in references [35,36]. It requires specification of grid point locations along 
the boundary - in this case, both the inner boundary (airfoil surface) and the outer 
boundary. The solution algorithm is outlined below. 

The Poisson equations are used to generate a boundary-fitted, curvilinear 2-D 
grid: 

(n £yy 

Vxx d” tyyy 

where (£,7?) represent coordinates in the computational domain, (x,j/) represent co- 
ordinates in the physical domain, and P and Q are source terms which control the 
grid point spacing in the mesh interior. The computational domain is rectangular 
and the grid points within it are evenly-spaced. To simplify the evaluation of the 
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derivatives and to ease the specification of the boundary conditions, the above equa- 
tions are transformed to, and solved in, the computational domain. To do this, the 
roles of the independent and dependent variables are interchanged, and the equations 
become 

ocX££ 2/3x^ v "j” ^x^ - — J ( Px ^ + Qxjy') 

- 2^2/^ + 72/™ = ~J 2 { P V(. + QVv) (3- 1 ) 

where 

a = x l + yl 

P = x i x r , + y i y n 
7 = x l + y 2 
J = ~ x nVi 

All derivatives are approximated using standard second-order accurate finite dif- 
ferences. The spatial increments A£ and Ai 7 can, without any loss in generalization, 
be assumed to be constant everywhere and equal to 1. The grid point locations on 
the boundary must be specified, and an initial guess for the interior grid values must 
be made. Any number of standard relaxation schemes can be used to solve the sys- 
tem of two elliptic equations. Here, a fully implicit scheme - an alternating direction 
implicit (ADI) scheme - constructed using approximate factorization to convert the 
solution process to two one- dimensional inversions is employed. The solution method 
is discussed in more detail in reference [37]. 

This method with P = Q = 0 provides no control over the grid point spacing near 
a boundary. The grid points tend to be pulled away from the surface by the Laplacian 
elliptic solver. Sorenson and Steger [38] developed a technique for defining P and Q 
such that the angle at which the £ = constant grid lines intersect the boundary, as 
well as the distance between the boundary and the first off-boundary grid point on 
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these grid lines, could be specified. In this way, grids having a very fine grid spacing 
near surfaces (for viscous calculations) can be generated. Also, orthogonality of the 
grid at the boundaries can be specified, if desired. 

In this current study, the elliptic grid solver has been used to generate the grids 
for the 2-D and 3-D freestream calculations. Since Sorenson and Steger’s boundary 
control technique has not as yet been implemented (i.e. P = Q = 0), post-processing 
in the form of exponential stretching was employed to obtain the desired grid density 
in the boundary layer region. 

For the wing/rotor interaction computations, a grid with a flat outer boundary 
in the plane of the rotor is desired so as to enable an easier and a more accurate im- 
plementation of the boundary conditions that correspond to the rotor. An outer grid 
is then defined to encircle the inner zone. See , for example, Figures 5.20 and 5.21. 
These grids have been generated purely algebraically. Elliptic grid smoothing has 
not as yet been applied to the rotor /wing meshes. 
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Chapter 4 

Boundary Conditions 


4.1 General Remarks 

The finite difference solution of the Navier- Stokes equations requires specification of 
boundary conditions on all domain boundaries. In the numerical method employed 
in this study (described in Chapter 2), the boundary conditions are applied explicitly, 
i.e. the flow variables at the boundaries are evaluated at the time level preceding that 
at which the implicit solution on the mesh interior is found. This permits greater 
flexibility in applying the boundary conditions to a variety of geometries and flow 
situations. At all grid points located on the mesh boundaries, each of the five flow 
parameters that make up the vector of conserved quantities, Q, must be updated 
explicitly - either by specifying them or by extrapolating from computed interior 
values. Referring then to the definition of Q in Section 2.2, the density the three 
components of mass flux /m, pv , and pw, and the total energy per unit volume, e, 
must all be updated at each time step. 

Determination of the boundary conditions representative of a lifting rotor require 
a separate analysis and will be discussed in a later section. First, those boundary 
conditions not pertaining to the rotor will be described. 
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4.2 Non-Rotor Boundary Conditions 


At grid points on the wing surface, for viscous (Navier- Stokes) calculations the no- slip 
boundary condition is imposed, i.e. all components of velocity are zero (it = v = 
w = 0). Inviscid computations can also be performed using the present computer 
code by omitting all viscous terms, i.e. the G v flux vector, and by applying the 
inviscid boundary condition on the wing - zero normal velocity. In the rectangular 
computational domain, this condition is easily satisfied by setting the contravariant 
velocity component normal to the surface, W , to zero. 

The pressures on the wing surface are found by solving the normal momentum 
equation (refer to [33,34]). The normal momentum equation is derived by taking the 
dot product of the vector comprised of the transformed x-, y-, and z- momentum 
equations, and the unit normal vector, n. The viscous effects are assumed to have an 
insignificant effect on the pressure at the surface and are neglected (typical boundary 
layer assumption). 

J x—mom i + y—mom j 4- z—mom k J ■ n = normal momentum 

where 

~ _ VC _ Cx ? + Cyj + (zk 

" ” |vcj - y/Q + q + a 

From the momentum equations, it can be seen that the normal-momentum equation, 
at the body surface reduces to dp/dn = p n = 0. Performing the above operations, 

Pi (CxCr + CyCy + Uz) + Pr, (*?xCx + VyCy + VzQ + Pi {(I + Cy + C) 

+ pU (Cx«€ + CyV{ + + pV (CxU n + (yV v + Cz^r,) 

= Pn \J Q + q + Q 

where U = V = 0 for viscous flow calculations, p ^ can be approximated by second- 
order one-sided differences and p ^ and p v by second-order central differences at the 
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surface. Re-arranging the equation, and applying approximate factorization, re- 
sults in an implicit solution algorithm for p at the surface which involves two one- 
dimensional tridiagonal inversions - one in the ^-direction and the other in the 77 - 
direction. Obtaining surface pressures using the above method yields a more accurate 
and stable solution method than simply using zero-order extrapolation. 

Once the surface pressures are known, assuming adiabatic conditions at the sur- 
face (i.e. no heat flux - zero normal temperature gradient), then the density is 
obtained by a zero-order extrapolation from the value at the nearest off-body point. 

The final quantity required, the total energy per unit volume, e, is calculated 
directly from Eq. 2.12 using all previously-defined quantities. 

Only one-half of the tilt rotor configuration is modeled due to the symmetry that 
exists in hover. Two parallel 2-D grid planes straddle the plane of symmetry and 
these grid points are used to specify the following symmetry boundary conditions 


Pi 

= P3 

(Hi 

= (H 3 

(Hi 

= - (H v 

(/Hi 

11 

ei 

= e 3 


This ensures that, at the centerline, there is no temperature gradient (assuming adi- 
abatic conditions) or pressure gradient, and that the spanwise component of velocity 
is zero. 

The 3-D Navier-Stokes code was modified to allow the computation of two- 
dimensional flows. Because fourth-order spatial differencing is used on the right 
hand side of Eq. 2.23, at least five parallel and identical 2 -D grids are required. For 
“pseudo 2-D” calculations, then, symmetry conditions are also applied between grid 
planes 3 and 5. This ensures that the code effectively sees an infinitely long wing 
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having constant airfoil section. 

On all outer boundaries other than the plane of symmetry, inflow or outflow 
boundary conditions are specified. The flow is essentially inviscid in these far- field 
regions. The Euler equations are hyperbolic partial differential equations. Applying 
a method of characteristics analysis to hyperbolic PDE’s helps to determine the 
appropriate boundary conditions for inflow and outflow boundaries. For subsonic 
flow in three dimensions, four of the characteristic velocities are positive and the 
fifth one is negative. For a subsonic inflow boundary, then, four independent flow 
thermodynamic and kinematic flow properties should be specified, and one should be 
extrapolated from the interior of the flow domain. For a subsonic outflow boundary, 
on the other hand, only one property should be specified and four extrapolated. 
For freestream calculations for flow at an angle of incidence of -90 degrees, inflow 
boundary conditions are applied at the grid points on the outer boundary that lie 
above the chord plane. Freestream Mach number and flow angle are specified. Also 
the pressure is set to be freestream ambient. The density is extrapolated from the 
interior using zero-order extrapolation. The outflow boundary is defined as those 
points on the outer boundary which lie below the chord plane. Here zero-order 
extrapolation from the nearest interior point is used for four flow properties. The 
pressure is specified and assumed to be again at freestream ambient. 

The grid points of the O-grid along the airfoil surface (j index) wrap around the 
airfoil and overlap by three points. This couples the j = 1 and the j = J M AX 
boundaries, which in the rectangular volume computational domain, are at oppo- 
site ends. The following periodic boundary conditions are then imposed, where the 
subscripts are for the j index: 


<h 

— <lJMAX-2 

Qj MAX 

— 93 


where q represents an element of the vector of conserved quantities, Q , and the above 
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is applied to all five elements. Also, for consistency, 

<72 = qjMAX-\ 


is specified. 

4.3 Modeling the Rotor 

4.3,1 Approach 

As discussed in Chapter 1, detailed modeling of individual rotor blades using CFD 
is a formidable task requiring computer resources that push the currently-available 
technology. Because the focus of this work is wing download, the problem is rendered 
more tractable by employing a simpler model for the rotor. The rotor is assumed 
to be an actuator disk. The combined position- and time- dependent effects of the 
blades on the flow are averaged. In the analysis presented here, radial variations in 
loading are permitted. On any given annular ring of the actuator disk, the load is 
assumed to be uniform. 

A classical analysis of propellers/ rotors is employed which combines axial and 
angular momentum conservation through the rotor disk with blade element strip 
theory. The latter allows for the incorporation of empirical data which can introduce 
viscous (Reynolds number) effects and Mach number effects at the rotor, if desired. 
The current analysis presented below does not include the effect of the wing on the 
rotor. Clark [8], using a panel method applied to the tilt rotor configuration, showed 
that although the local effects of the wing on the rotor are significant, the global 
effects on total torque and power are surprisingly small due to counterbalancing 
local effects. 

The analysis below is used to independently obtain the pressure rise through the 
rotor disk and the swirl velocity immediately downstream of the disk. These are 
then used as boundary conditions for the Navier-Stokes code. 


34 


4.3.2 Combined Momentum Conservation - Blade Element 
Analysis 

Glauert [39], McCormick [40], and Prouty [41] provide good discussions on momen- 
tum conservation and blade element analysis applied to propellers and/or rotors. 
The theory, as it is applied here, is outlined below. 

The analysis assumes that (i) the slipstream produced by the rotor does not 
contract. This is a valid assumption for lightly or moderately loaded rotors, which 
is assumed to be the case here, (ii) There is no radial interference, i.e. the flow is 
two-dimensional at every spanwise station of the rotor. This assumption is valid 
except for a very small region near the blade tip. Small empirical corrections can 
be made for this effect, if desired, (iii) The torque and axial loads produced by 
the blade elements of the finite number of blades at distance r from the rotation 
axis, are averaged over the entire annulus, (iv) The rotation axis is aligned with the 
freestream u 0 0 flow. We assume here, for convenience, that the propeller /rotor axis 
is horizontal, (v) The flow is incompressible. 

Expression for the Thrust 

Applying momentum conservation to an annulus of a cylindrical control volume at 
radius r, having width dr , and which extends to infinity upstream and downstream 
of the rotor, 

dT = p (27 rr dr) (Uqq + u> 0 ) (2i«a) (4.2) 

mass flow velocity change 

where dT is the average elemental thrust force in the axial direction on an annulus 
of area 27r rdr, u <*, is the freestream velocity, w a is the induced axial velocity in the 
rotor disk, and 2 w a is the induced axial velocity far downstream. Eq. 4.2 equates 
the net axial force with the rate of change of momentum (the mass flow multiplied 
by the net change in axial velocity). 

That the induced axial velocity far downstream is twice that at the rotor plane 
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Figure 4.1: Relative velocity and force diagram at a rotor blade element. 

can be shown in two ways, (i) Simple actuator disk theory using the principles of 
mass continuity and momentum conservation, and applying the Bernoulli equation 
twice — once upstream of the rotor disk, and once downstream of the rotor disk, (ii) 
Model the circulation on the rotor blade and the vorticity shed into the wake by 
horseshoe vorticies placed at each blade element. 

The elemental thrust force can also be expressed in terms of the 2-D lift and 
drag of the airfoil at a blade element. Referring to the relative velocity and the force 
diagram of Figure 4.1, it can be seen that 

dT = B ( dL cos {(f> + a,) — dD sin ( <f> + a,) ) (4-3) 


where the total element thrust dT is equal to the sum of the contributions of all 
blades, B is the total number of blades, a, is the induced angle of incidence, and <f> 
is defined in Figure 4.1. 

From the definition of the aerodynamic force coefficients, 


dL 


2 P ^ C< ^ r 
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dD 


( 4 . 4 ) 


= 2 P u ln C d c dr 

where C/ and Cj are the two-dimensional airfoil lift and drag coefficients, respectively, 
c is the local blade chord, and u e jj is the effective velocity that the blade element 
sees. 

Substituting Eqs. 4.4 into Eq. 4.3 , 

dT = ]- pu 2 e jj B cdr (Ci cos ( <f> -f a,) — C d sin (<j) + a,) ) (4.5) 

Equating Eqs. 4.2 and 4.5, and observing from Figure 4.1 that 

(^ + w a ) 2 = u\ J} sin 2 (4> + a;) 

we get: 

W a _ Be (Cl cos (<ft -f a,-) - Cd sin (<ft -f ap) ^ ^ 

Uoo + w a S-Jtr sin 2 (<f> + Oj) 

Note that the aerodynamic coefficients can be expressed as a function of angle of 
incidence, a, i.e. 


C, = /i(or) 

= h(a) 

where, from Figure 4.1, 

q = 0 + 6 - ((f) + ai) (4.7) 

where j3 is the blade pitch angle setting, and 6 is the local blade twist angle relative 
to the twist at the 75% radial station. 

Empirically-derived section lift and drag characteristics as a function of angle 
of incidence can therefore be used in this analysis, thereby improving the predicted 
results for the rotor, and at the same time reducing the computational effort that 
would otherwise be required to generate this data using a Navier-Stokes solver, for 
example. 
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Note that for hover, u <*, = 0, and the left hand side of Eq. 4.6 becomes unity. 
The only unknown in this equation, then, is a,. The angle 4> is zero when u Q 0 = 0 
(see Figure 4.1). 

Expression for the Torque 

The swirl velocity, w t , is zero ahead (upstream) of the rotor, w t at the rotor, and a 
constant 2w t from immediately downstream of the rotor face to infinity downstream, 
in the inviscid analysis employed here. 

The conservation of angular momentum equates the net torque produced by an 
annular element of the rotor on the fluid, with the rate of change of angular momen- 
tum (the mass flow multiplied by the net circumferential change in velocity multiplied 
by the moment arm r). 

dQ = p (2nr dr) (uqq + te a ) (2u^) r (4.8) 

mass flow velocity change 

Similar to the thrust analysis, the elemental torque contribution averaged on an 
annulus of width dr located a distance r from the axis of rotation can be written in 
terms of the elemental lift and drag loads. Referring to Figure 4.1, 

dQ = B ( dL sin (<f> + a,) + dD cos (<j> + aj ) r (4-9) 

Substituting Eqs. 4.4 into Eq. 4.9 and equating to Eq. 4.8 yields, upon re-arrangement, 

wt _ Be (Ci sin (<j> + aj + Cd cos (<£ + Qp ) ^ ^ 

Uoo + w a Snr sin 2 ( <j> + a,) 

In hover, when u <*, = 0, the left hand side of Eq. 4.10 reduces to w t /w a . Since a, is 
determined from Eq. 4.6, then Eq. 4.10 can be solved for w t /w a . 

Calculation of Required Quantities 

Referring to Figure 4.1, it can be seen that 

W = yj (Uoo + Wof + K - W<) 2 sin a * = \J W l + w t 
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where u> is the angular frequency of rotation of the rotor. Re-arranging and setting 
= 0 for hover, we get for the local induced axial velocity 


ur 


W n — 


(4.11) 


^cos^. + K/^ ( w J Wa ) 
sm or, v w 

Once w a is calculated using Eq. 4.11 knowing w t /w a from Eq. 4.10, then the swirl 
velocity w t is immediately known as well. 

Applying the Bernoulli equation upstream of the rotor and again downstream of 
the rotor, we find that the pressure rise across the rotor, A p, in hover, is 

Ap = 2 pw 2 a 

The total thrust and torque on the rotor disk is calculated by integrating the radial 
thrust gradient and radial torque gradient, respectively. Eqs. 4.2 and 4.8 can be 
easily written in terms of dT/dr and dQjdi\ respectively. Then 


Tjot = f (dT/dr) dr 
Jo 

Qtot = f (dQjdr)dr 
Jo 


The total thrust coefficient and total torque coefficient are calculated using the stan- 
dard definitions: 

Ttot 


Cj — 
C Q = 


P A { U lo + ( w R ) 2 ) 

Qtot 


P A + (^^) 2 ) R 
where A is the rotor disk area and R is the blade radius. 


Solution Procedure 

For a given rotor having specified blade geometry and for a given thrust coefficient, 
Cti an iterative solution procedure outlined below gives the radial distributions of 
axial and swirl velocities and pressure rise, at a converged value of blade pitch angle 
setting, /?. 
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1. Input rotor RPM, and rotor geometry including: 

• number of blades 

• rotor radius 

• blade planform shape, i.e. local chord, c, distribution 

• blade twist relative to the 75% radial location 

2. Input 2-D airfoil characteristics: C/, Cd vs. a. 

3. Input desired total thrust coefficient, Cx- 

4. Make an initial guess for the blade pitch angle, /?. 

5. Divide the rotor disk into a number of annular rings. At each, calculate the 
induced angle of incidence, o,, using Eq. 4.6. Calculate the effective angle of 
incidence, a, from Eq. 4.7. Use the inputted 2-D airfoil data to determine the 
Ci and Cd corresponding to a. 

6. Calculate w a , w t , A p, dT/dr, and dQ/dr at each annulus. 

7. Integrate along the radius to obtain Cx and Cq. 

8. If Cx is within a certain specified tolerance from the desired, inputted value of 
Cx , then the solution is converged; otherwise, go to (4) and repeat steps (4) 
through (8). 

Figure 4.2 is an example of the output generated from this analysis, for a 3- 
bladed, 7-ft diameter rotor similar to the 0.16 scale model tested at the NASA Ames 
Outdoor Aerodynamic Research Facility. For this example, NACA 0012 airfoil data 
was substituted for the aerodynamic characteristics of the actual rotor blade sections. 


A 
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Figure 4.2: Sample output from the momentum theory/blade element analysis for a 
7-ft diameter, 3-bladed rotor at 2090 RPM and with Ct = 0.0065. 

4.3.3 Grid Used for Modeling the Rotor 

To limit the size of the computational domain, it was originally intended to extend 
the inflow boundary only to the horizontal plane containing the rotor. Solution 
convergence problems were encountered, however, due to the close proximity of the 
boundary to the wing whose presence induces large flow gradients near the inflow 
plane. It was decided, therefore, to increase the size of the physical domain so 
as to include the effect of the large-scale entrainment of flow back into the top of 
the rotor. An inner zone surrounding the wing and extending to the rotor plane 
was developed. An outer zone was also created which encircled the inner grid and 
extended approximately 20 rotor diameters from the wing in all directions. The 
outer grid line of the inner zone matches exactly the inner grid line of the outer 
zone. See, for example, Figures 5.20 and 5.21. Consecutive, independent solutions 
are performed on the inner and outer zones until adequate conveigence is achieved. 

Everywhere along this common boundary between the two zones except foi points 
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located in the region defined by the rotor disk, simple averaging of the flow properties 
is performed to couple the two solutions. For a given spanwise grid station, k, and 
a given chordwise index, j , the most recently computed values of the flow properties 
at the interior grid points immediately to either side of the common boundary are 
averaged. This same procedure is applied both for the outer boundary points of the 
inner grid and the inner boundary points of the outer grid. This averaging assures a 
continuity of mass across the common boundary. 

The outer boundary grid points of the inner zone that lie within the rotor disk are 
viewed as an inflow boundary. From a method of characteristics analysis, as discussed 
in Section 4.2, a subsonic inflow boundary requires specification of four quantities and 
extrapolation of the fifth. The interior values of pu, pv, and pw are averaged across 
the disk as described above, and these averages used as inflow boundary conditions 
for the inner zone. The pressure, p, is obtained from the pressure at the first interior 
point of the outer zone (just above the rotor face). To this is added the pressure 
rise, A p, obtained in the momentum conservation theory/blade element analysis of 
Section 4.3.2 . The density, p, is extrapolated from the interior of the inner zone 
using zero-order extrapolation. To include the effects of swirl, the tangential velocity, 
u>t, obtained from the analysis of Section 4.3.2, can be decomposed into a u and v 
contribution in the horizontal plane and used to replace the values above obtained 
by averaging. 

For the inner boundary grid points of the outer zone that lie within the rotor 
disk, search for the best combination of boundary conditions is still being pursued. 
Currently, averaging across the rotor is used to specify the boundary conditions for 
all the flow variables but the pressure. For the pressure, the prescribed pressure rise, 
A p, is subtracted from the most recently-computed value of p at the nearest interior 
point in the inner zone. Specifying A p across a zonal boundary was first suggested 
to the authors by Tavella [42]. 
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Figure 4.3: One-point overlap grid for wall jet boundary conditions. 

4.4 Wall Jet 

In previous CFD studies (for example, reference [19]), multiple zones were used to 
model a tangential, circulation control jet on the surface of an airfoil. Preliminary 
calculations during the course of the current research, however, have demonstrated 
that a single-zone grid, an approach suggested by Tavella [42], is suitable for this 
computation. 

The grid must have sufficient grid density in the jet region to resolve the jet flow 
and associated entrainment of the outer flow. For easier implementation of boundary 
conditions for the wall jet, the 3-point grid overlap in the chordwise j index normally 
used, is reduced to a 1-point overlap (see Figure 4.3). 

The wall jet boundary conditions are specified on the j — JMAX line. The grid 
is defined such that 10 - 20 grid points lie within the ~ 0.1% chord jet slot height. 
Beyond this jet region, the values of the flow properties at j = 1 and j = JM AX are 
obtained by taking the average of the values at the corresponding interior points at 
j — 2 and j = JMAX—1. Within the slot region, the j = JMAX grid points are viewed 
as an inflow boundary. A total pressure corresponding to the plenum pressure used 
in the model tests, a total temperature assumed to be freestream ambient, and jet 
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exit angle are all specified. The jet is assumed to be tangent to the airfoil surface at 
the jet exit. The magnitude of the jet exit velocity is calculated assuming isentropic 
expansion of the compressed air from the plenum to local static pressure at the jet 
exit . The j = 1 grid points that overlap the j = JM AX points in the jet are viewed 
as outflow boundary points. Pressure is specified and set equal to that at the jet 
exit. Zero-order extrapolation from the j = JMAX — 1 values is used for all the other 
flow properties. 

4.5 Initial Conditions 

The steady-state (or pseudo steady-state in this case where some residual unsteadi- 
ness due to vortex shedding off the wing leading and trailing edge may exist) solution 
of a hyperbolic system of partial differential equations is independent of the initial 
conditions. In our case, the interior points can be seeded initially with a freestream 
flow close in magnitude to the expected rotor slipstream velocity, or the interior could 
be set initially to have zero flow everywhere. Both of the above initial conditions 
yield the same final results, with similar rates of convergence. 
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Chapter 5 

Discussion of Results 


5.1 Preliminary Comments 

The computation of tilt rotor flowfields using the Navier-Stokes equations is a chal- 
lenging task. A problem such as this is best divided into smaller, more manageable 
parts, and a step-by-step approach used to reach the final objective. That is why 
this preliminary CFD study of tilt rotor flowfields has focused on the calculation of 
subsets of this complicated 3-D flow. 

In the following discussion, freestream results both in 2-D and in 3-D are presented 
for an airfoil and a wing, respectively, at an angle of incidence of -90 degrees. The 
rotor alone was modeled as an actuator “line” in 2-D and as an actuator disk in 3-D 
with and without swirl. These results are followed by two-dimensional wing/rotor 
interaction computations. Preliminary results for 3-D wing/rotor interaction are also 
discussed. Finally, the results of a 2-D tangential blowing calculation are presented. 
These computations serve as building blocks necessary to ensure the successful com- 
pletion of the final objective - to compute accurately the three-dimensional tilt rotor 
flowfield with a wing having tangential blowing for leading edge separation control. 

All plots shown in this chapter were generated using PLOT3D, an interactive 
graphics package developed at NASA Ames Research Center [43,44], 
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5.2 Freestream Calculations (No Rotor) 

5.2.1 Results in Two Dimensions 

Figure 5.1 shows two views of a typical grid used to calculate the freestream flow 
around the V-22 airfoil section at an angle of incidence of -90 degrees. In order to 
resolve the fine details of the flow in the region of a blunt trailing edge, very high 
grid resolution is required. The focus of this work was not on the trailing edge flow 
alone. For this reason, as is often done in CFD airfoil work, the original trailing 
edge is extended to a point and the modified airfoil is re-scaled to the original chord 
length. The V-22 airfoil trailing edge thickness is considerably less than 0.5% chord 
so the modification to the airfoil section is slight. A 2-D elliptic grid solver is used to 
obtain a smooth distribution of grid points. Exponential clustering produces a fine 
mesh near the airfoil surface. As discussed in Section 4.2, five identical, parallel 2-D 
grid planes are required for the two-dimensional calculations so as to accommodate 
the 5-point, fourth-order differencing stencil used on the right hand side of Eq. 2.23. 
Symmetry boundary conditions are imposed on both ends of the grid, to simulate 
an infinitely-high aspect ratio wing. 

Figures 5.2 to 5.5 are all results of a “pseudo” 2-D calculation. A freestream flow 
of M = 0.2 at a = —90° is imposed on the upper inflow boundary. A Reynolds 
number of 0.5 x 10 6 is selected as being representative of the small scale (~ 0.16 
scale) model flow in tests that were - and that will be - undertaken at NASA Ames’ 
Outdoor Aerodynamic Research Facility (OARF) - refer to Section I.2.T. A time- 
accurate computation is performed, i.e. the time step is sufficiently small to resolve 
the unsteadiness in the flow. 

Figure 5.2 shows the velocity vectors around the airfoil at two different points 
in time. In the upper picture, a vortex has just been shed off the airfoil leading 
edge. In the lower picture, a vortex has just been shed off the trailing edge. This 
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(b) Close-up view of the airfoil gridding. 



Figure 5.1: Views of a 2-D O-grid around the V-22 airfoil. 
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asymmetric vortex shedding is typical of the vortex streets created behind bluff 

bodies at these Reynolds numbers. The frequency of this shedding was computed to 

be about 16 cycles per second. For comparison purposes, a 2-D cylinder, at the same 

Reynolds number, whose diameter is equal to the airfoil chord, typically produces 

vortex shedding at the rate of about 30 cycles per second. This value was obtained 

from an empirical relation in reference [45]: 

fd 19.7 

— = 0.198 1 - — - 
Uoo \ He 

where / is the frequency of vortex shedding and d is the cylinder diameter which, 
for the above calculation, was set to 0.447 meter which is the 0.16 scale model wing 
chord. No experimental data exists for the V-22 airfoil to confirm the accuracy of 
this unsteady, computational result. Examining the velocity profile at the airfoil 
surface, it can be observed that the boundary layer has been resolved. Figure 5.3 
is a close-up view of the flow around the wing leading edge. Separation is seen to 
occur at a location on the airfoil upper surface slightly aft of the leading edge. This, 
of course, changes somewhat with time due to the unsteady nature of this 2-D flow. 
Note that in this and all subsequent calculations discussed in this report, laminar 
flow was assumed - i.e. no model has as yet been implemented to model turbulence 
effects . 

Figure 5.4 is a Mach contour plot of the flowfield at one point in time. Note that 
the flow stagnates in the mid-chord region on the upper surface, and beyond the 
leading and trailing edges, the flow accelerates to about M = 0.26. 

A computation was performed where large time steps were taken to accelerate 
convergence to a “pseudo” steady-state solution. The resulting pressure distribution 
is shown in Figure 5.5. The general shape of the C p distribution compares favorably 
with the download distribution of Figure 1.3. The loop in the distributions near the 
leading edge indicates a deceleration of the flow even before it gets to the leading edge. 
This was observed previously in Figure 5.3. The base pressure is fairly uniform along 
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Figure 5.5: Surface pressure coefficient distribution on the V-22 airfoil, in a 

A/*, = 0.2, a = —90° flow. 


the airfoil lower surface; this is also seen experimentally. The computed base pressure 
is higher than wind tunnel measurements. This is most likely due to the inaccuracy 
of the steady-state calculation for time dependent problems. The computed flowfield 
below the airfoil showed two symmetric vortices of equal strength but opposite sense 
of rotation whose centers were located over a chord length downstream of the airfoil. 
Certainly, in the future, it would be better to run the code in a time accurate fashion 
and take a time average of the desired flow characteristics. A pressure integration 
around the airfoil yields a computed 2-D drag coefficient Cd of 1.2. This compares 
to Cd ~ 1.6 obtained from wind tunnel tests of the XV- 15 airfoil [6]. 

5.2.2 Results in Three Dimensions 


As discussed in Section 3.1, a grid for three-dimensional calculations is generated by 
placing a series of 2-D grids at various spanwise stations along the wing and beyond 
the wing tip. Exponential stretching increases the grid density near the tip where the 
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three-dimensional effects are most significant. The wing is based on the dimensions 
of the 0.16 scale model used in the NASA Ames OARF tests. It has an aspect ratio 
of 6.1 and is of constant chord with no twist or sweep. By comparison, the V-22 wing 
has an aspect ratio of 5.5 and 6° of forward sweep. The outer boundary of the grid 
extends 20 chord lengths in all directions. As mentioned previously, due to symmetry 
of the flow about the centerline, only one half of the wing needs to be modeled. For 
the results shown in this section, the grid has 47 points defining the wing surface. 
Each grid line extending from the body to the outer boundary is defined using 33 
points, with exponential clustering near the surface. There are 35 parallel 2-D grid 
planes defining the spanwise grid distribution. This yields just over 54,000 points 
in total. A converged solution takes about 700 time steps and about 30 minutes of 
CPU time on the Cray-XMP supercomputer. Figure 5.6 shows two views of the grid 
employed. Figure 5.6(a) is a cutaway perspective view of the grid, and Figure 5.6(b) 
is a horizontal cut through the wing chord (x — y ) plane and including the upper 
surface gridding. The wing tip is defined arbitrarily to be the last 5% of the wing 
semi-span, and it has a chord variation that is elliptic, with each tip airfoil section 
having the same thickness-to-chord ratio, tjc. It can be seen from Figure 5.6(b) that 
the tip region is not well-defined, and the mesh is highly skewed here. Improvements 
in gridding the tip region will be pursued in future work. 

Figures 5.7 and 5.8 are both views of the flow features in a vertical plane running 
spanwise through the wing mid-chord. Figure 5.7 is a plot of velocity vectors for a 
freestream flow of M = 0.2 at a = —90°, the same as for the 2-D case. Figure 5.8 
shows the corresponding Mach number contours. Both of these plots show that 
much of the span of the wing in the mid-chord region sees stagnated flow. The flow 
becomes increasingly more spanwise as the tip is approached, as expected. The flow 
accelerates around the tip, and viscous effects give rise to a very pronounced vortex 
beneath the tip. 
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(a) Cutaway view of the grid. 



x/c 


(b) An x — y cut through the wing chord plane and including the wing upper 

surface grid distribution. 


Figure 5.6: Views of a 3-D grid around the finite-span wing. 
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Figure 5.9 is a view of the velocity vectors above the wing surface and in the x—y 
plane beyond the wing. The stagnation region runs from the centerline to about half 
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Figure 5.9: Velocity vectors just above the wing and beyond the wing in the x—y 
plane, in a — 0.2, a = —90° flow. 

way to the wing tip. In this region, the flow is primarily two-dimensional in character. 
From the mid-span region outboard to the tip, the flow becomes increasingly spanwise 
before rounding the tip where it is rapidly re-directed by the freestream. 

Figure 5.10 shows five separate chordwise surface pressure coefficient distribu- 
tions, each at a different spanwise location on the wing. Three-dimensionality is 
again observed in the form of changing pressure distributions particularly in the tip 
region. For all spanwise stations, the wing leading edge induces a greater flow ac- 
celeration than the trailing edge - as evidenced in lower peak C p 's at the leading 
edge. 

Figure 5.11 shows seven separate upper surface pressure coefficient distributions, 
each for a different percent chord location. Here, tip effects are clearly shown. The 
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Figure 5.10: Typical chordwise surface pressure coefficient distributions at various 
spanwise stations along the wing, in a M 0 0 = 0.2, a = —90° flow. 

pressure coefficient for all but the leading edge C p distribution drops sharply at the 
tip. Kinks in the tip region are due to insufficient grid- point density. Along the 
wing leading edge, as the tip is approached from the mid semi-span, the flow is 
accelerated not only by the 2-D flow effect around the leading edge, but also by 
the increasing, 3-D spanwise flow. That is why the negative pressure coefficient 
starts to increase well inboard of the tip. After a certain point, as evidenced by the 
peak in the leading edge C p distribution, the spanwise flow dominates and the flow 
acceleration diminishes because less of the flow travels around the leading edge to 
the lower surface. 

5.3 Rotor Alone in Two and Three Dimensions 

Initially the inflow boundary (upper boundary - see, for example, Figure 5.6) for 
the wing/rotor computation was to be the plane of the rotor. The rotor was to 
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Figure 5.11: Spanwise upper surface pressure coefficient distributions at seven dif- 
ferent chordwise locations along the wing, in a A/^ = 0.2, a = —90° flow. 
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be modeled by specifying a pressure and velocities corresponding to a lifting rotor 
(obtained from the momentum theory /blade element analysis of Section 4.3). It is 
unclear, however, what boundary conditions to impose at those grid points that are 
on the inflow boundary but away from the rotor disk. Being an inflow boundary, 
pressure and velocity should be specified, for consistency. These, however, are not 
known to sufficient accuracy. A previously- developed vortex panel model of the tilt 
rotor was deemed inadequate to provide this information. Additionally, no velocity 
or pressure surveys have been performed that could be of assistance. Computational 
results that were obtained by extrapolating values of the flow properties from the 
interior of the solution to the upper boundary were not satisfactory. Another draw- 
back of this initial approach is the fact that because the pressures and velocities were 
specified at the rotor face, the influence of the wing was completely decoupled from 
the rotor. 

For these reasons, it was decided to model the rotor using two zones, and a 
multiple zone algorithm. The common boundary between the two zones is selected 
to coincide with the location of the rotor. This allows specification of a pressure rise, 
A p , across the disk without having to fix the static pressure itself. The pressure 
in the rotor plane, then, can float under the influence of the wing below. The use 
of two zones also extends the solution domain much further thereby permitting the 
resolution of the fountain flow and the large scale recirculation pattern through the 
rotor. 

Figures 5.12 and 5.13 show the results of two different two-dimensional calcu- 
lations using a Cartesian grid comprised of two zones, one above the other. Each 
zone has 47 grid lines in the direction parallel to the rotor axis of rotation and 17 
grid lines perpendicular to the rotor axis. The common boundary between the zones 
is used to define the rotor boundary conditions. Exponential stretching is used to 
cluster the grid points in the location of the rotor. The freestream velocity was set to 
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M 0 o = 0.001 so the flow is driven exclusively by the specified pressure rise across the 
rotor. Here, the magnitude of the pressure rise, A p, is 10% of freestream ambient 
pressure. The actuator “line” of these 2-D calculations is assumed to be uniformly 
loaded, i.e. A p = constant. 

Figure 5.12 shows the instantaneous particle traces (same as streamlines in steady 
flow) through the rotor disk. For inviscid flow, the flow accelerates through the 



Figure 5.12: Particle traces released in a flow through a 2-D actuator disk, where 
A p = O.Olpoo. 

rotor disk reaching a maximum velocity downstream of the disk. Theoretically, this 
flow momentum is maintained indefinitely. The slipstream, therefore, contracts to 
a minimum width below the rotor, and thereafter remains unchanged. Examining 
Figure 5.12, it is seen that the computed slipstream contracts to a minimum width, as 
expected, but then begins to expand slowly indicating that there is a corresponding 
deceleration of the flow. This behavior is unrelated to the contribution of the viscous 
terms in the Navier- Stokes equations because the same results are obtained for the 


59 



Euler equations alone. Admittedly, an actual rotor slipstream diffuses within a few 
rotor diameters. Since it is not the physics modeled by the equations that is causing 
the diffusion, it must, therefore, be due entirely to numerical dissipation. Numerical 
dissipation is caused by the truncation error associated with the spatial differencing 
of the partial differential equations. It increases with grid cell size. Mesh refinement, 
then, i.e. a finer mesh, would reduce the computed diffusion of the rotor slipstream. 

Figure 5.13 shows the pressure above and below the rotor. As anticipated from 



Figure 5.13: Pressure contours for a flow through a 2-D actuator disk, where 

A p = O.Olpco. 

simple momentum theory, the pressure far above the rotor is at freestream ambient 
(non-dimensionalized = I/7 = 1/1.4 = 0.7143). The static pressure drops to a 
minimum value just above the rotor due to the accelerating flow being drawn into 
the rotor disk. The rotor produces an increase in total pressure, which in this figure, 
is seen as a step in the static pressure immediately downstream of the rotor. The 
pressure then decays back towards its original freestream value as the flow accelerates 
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to a near constant value downstream of the rotor. As expected, the flow gradients 
increase as the rotor disk is approached. 

Three-dimensional results are obtained for a rotor by creating a mesh in the form 
of a rectangular box. A horizontal plane through the middle of this grid forms the 
common boundary between the upper and lower zones. Each zone is comprised of 
35 stacked, parallel 2-D grid planes. Each grid plane is identical to that used for the 
2-D calculations discussed above. 

Figures 5.14 - 5.19 are results from a computation where A p across the rotor 
is set to a uniform O.OSpoc, the freestream is essentially zero ( M <*, = 0.001), and 
the rotor diameter is 6 units. The plots focus on a small portion of the entire 
domain which extends to 100 units in all directions. Figure 5.14 shows the computed 
velocity field in the near vicinity of the lifting rotor on a vertical plane through 
the rotor axis. Note that far above the rotor, the flow velocity is near zero. The 



Figure 5.14: Velocity vectors in a vertical plane through the 3-D actuator disk, where 
A p = 0.05poo. 
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flow accelerates continuously down through the rotor until it reaches a maximum 
value within a rotor radius downstream of the rotor face. Simple, incompressible 
momentum considerations in 1-D indicate that the maximum axial velocity, induced 
by an actuator disk in zero freestream flow, is twice that of the velocity in the 
rotor plane. Qualitatively, this can be seen in the figure by comparing the velocity 
vector lengths in the rotor plane with those far downstream. Figure 5.15 shows the 
accompanying particle traces through the rotor disk. Figure 5.16 is a plot of the 



Figure 5.15: Particle traces in a vertical plane through the 3-D actuator disk, where 
A p = 0.05poo. 

Mach number contours in a vertical plane through the center of the actuator disk. 
Note that the flow accelerates smoothly through the disk and the flow gradients 
increase in the vicinity of the rotor, as expected. With the given pressure rise of 5% 
of freestream ambient pressure, a simple analysis assuming incompressible, 1-D flow 
and applying the Bernoulli equation above and below the rotor yields a maximum 
induced Mach number of 0.267. This quick check provides confidence in the CFD 
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Figure 5.16: Mach number contours in a vertical plane through the 3-D actuator 
disk, where A p = 0.05poo. 

solution. As observed in the 2-D case, here also the rotor wake diffuses due to 
numerical dissipation. 

Figure 5.17 shows the computed Mach contours in the plane of the rotor. The 
squared contours in the outer regions are due to a combination of insufficient grid 
point density and the contour plotting algorithm itself. Certainly, to study a rotor 
alone, a cylindrical mesh is more appropriate. A rectangular type of mesh, however, 
fits more easily in with the rest of the tilt rotor configuration. 

Figures 5.18 and 5.19 show the effect of swirl on the flowfield. The tangential 
velocity (swirl) induced by the rotor is estimated from either experimental data or 
the momentum theory /blade element analysis discussed in Section 4.3 . For this 
test case, an arbitrary swirl Mach number of M = 0.05 is imposed uniformly on the 
upper boundary of the lower zone at the rotor location. The resulting flow pattern 
immediately downstream of the rotor is shown in Figure 5.18. It can be seen that the 
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Figure 5.17; Mach number contours in the rotor plane, where A p = 0.05p 
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Figure 5.18: Velocity vectors in a horizontal plane immediately downstream of the 
rotor plane showing the effect of swirl, where A p = 0.05poo and M swirl = 0.05. 










flow beyond the rotor is influenced by the specified rotational motion. It is obvious 
from this plot that a greater number of grid points is required to better represent 
the physical extent of the rotor disk. Figure 5.19 is a projection of the particle traces 
onto a vertical plane running through the center of the rotor and wake. Clearly, the 
effect of the swirl imposed at the rotor face is communicated to the rotor wake. 



Figure 5.19: Particle traces projected on a vertical plane through the actuator disk 
showing the effect of swirl, where A p = 0.05poo and M sw i r i = 0.05. 


5.4 Wing/Rotor Interaction 

5.4.1 Results in Two Dimensions 

With the experience gained from studying the wing and airfoil in a freestream, and 
the rotor alone modeled using two zones, the computation of 2-D airfoil/rotor inter- 
action is pursued. A more complicated grid is required. An inner zone is generated 
which extends from the airfoil surface to a flat boundary sufficiently wide to ac- 
commodate the complete rotor. Encircling this inner zone is another O-grid whose 
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inner boundary matches identically the outer boundary of the inner zone. The outer 
boundary of the outer grid extends about 20 rotor diameters in all directions - see 
Figures 5.20 and 5.21. The common boundary between the two zones is high- 



x/c 


Figure 5.20: Farfield view of the 2-D, two-zone grid used for computations of 
wing/rotor interaction. 

lighted. As in the freestream computations, 47 points are used to define the airfoil 
surface. In the inner zone, 22 points define each grid line normal to the surface. In 
the outer zone, 12 points define each grid fine. Again, for “pseudo 2-D” calculations, 
five identical, parallel grid planes are required for the computation. The current grid 
is produced algebraically, resulting in some undesirable kinks in the grid interior. 
Smoother grids, using the elliptic grid solver, await future grid refinement. 

Figures 5.22 and 5.23 show farfield and nearfield views, respectively, of the com- 
puted flowfield. The rotor diameter is 4.7S airfoil chords (determined from 0.16 scale 
model dimensions), and the uniform pressure rise across the rotor is A p = 0.05poc- 
The rotor diameter and the rotor height above the airfoil, when non-dimensionalized 
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Figure 5.21: Nearfield view of the 2-D, two-zone grid used for computations of 
wing/rotor interaction. 

by airfoil chord, all match the small scale model geometry of the past and projected 
future tests at NASA Ames Outdoor Aerodynamic Research Facility (OARF). For 
this hover calculation, M <*, was set to 0.001 . The flow is driven purely by the 
0.05poo that is imposed across the zonal boundary at all points that lie within the 
rotor disk. Note that the near-stationary flow far above the rotor is pulled down into 
the rotor disk. The flow impinges on the airfoil where it stagnates at the mid-chord, 
bifurcates, and accelerates around the leading and trailing edges. An asymmetric 
separated flow region is produced beneath the airfoil very similar to that seen in the 
freestream calculations (refer to Figure 5.2). Again, as before, the Reynolds number 
is 0.5 x 10 6 and although turbulent flow is expected, only laminar calculations have 
been performed. 

A moderately-loaded tilt rotor in hover produces an average pressure rise of 
around O.Olpoo across the rotor disk. Although good computational results have 
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Figure 5.23: Nearfield view of the velocity vectors for 2-D wing/rotor interaction 
where A p = 0.05poo and the rotor diameter is 4.78. 








been obtained with this magnitude of pressure rise (refer to subsequent discussions), 
a higher A p induces a higher Mach number flow which results in better solution 
convergence behavior for the compressible flow code. For program development pur- 
poses, then, a somewhat higher A p is often used. 

Figure 5.24 is a plot of the Mach number contours around the rotor and airfoil. 
The slight discontinuity in Mach contours through the rotor disk may be due to 



Figure 5.24: Mach number contours for 2-D wing/rotor interaction, where 

Ap = O.OSpoo and the rotor diameter is 4.78. 

an insufficiently converged solution, or insufficient grid point density in this region. 
Figure 5.25 clearly shows the drop in static pressure as the rotor is approached from 
above and the pressure jump across the rotor. The influence of the wing on the 
pressure distribution at the rotor is made more obvious by comparing this plot with 
Figure 5.13 for a rotor alone (where the pressure is uniform along the width of the 
disk). 

Figure 5.26 shows the static pressure distribution on the surface of the airfoil in 
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Figure 5.25: Pressure contours for 2-D wing/rotor interaction, where A p = 0.05poo 
and the rotor diameter is 4.78. 

the presence of the lifting rotor. Note that for the non-dimensionalization used in the 
flow solver, p <*, = 1/7 = 0.7143. The shape of the pressure distribution is very similar 
to that for an airfoil in —90° freestream flow (Figure 5.5). These results are from 
steady-state (i.e. large time step) computations of the flow code. A time accurate 
calculation, on the other hand, yields a base pressure distribution that changes as a 
function of time (depending on the relative positions of the leading and trailing edge 
vortices with respect to the lower surface). 

5.4.2 Results in Three Dimensions 

A 2-zone grid suitable for computation of three-dimensional wing/rotor interaction 
is generated by stacking many 2-D grids in parallel. Various views of the 3-D mesh 
are shown in Figures 5.27 to 5.30. Figure 5.27 is a farfield view of the grid showing 
the outer dimensions. Figure 5.2S shows more details including the wing surface, 
and the inner and outer zones at the wing centerline (y = 0). The actual location of 
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Figure 5.26: Pressure distribution on the surface of the airfoil for 2-D wing/rotor 
interaction, where A p = 0.05poo and the rotor diameter is 4.78. 



Figure 5.27: Farfield, cutaway view of the 3-D, two-zone grid used for computations 
of wing/rotor interaction. 
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Figure 5.28: Nearfield, cutaway view of the 3-D, two- zone grid used for computations 
of wing/rotor interaction. 

the rotor is superimposed on this view. Figure 5.29 is a horizontal cut in the plane 
of the rotor. The physical extent of the rotor disk is outlined here as well. Rotor 
boundary conditions are imposed at all points which lie on this plane and are within 
the rotor diameter. It is obvious from this figure that for accurate quantitative 
results, further grid refinement is required to define the rotor, particularly in the 
region outboard of the wing tip. Also, if the fountain effect is to be predicted, grid 
density must be increased in the centerline region of the mesh. The results presented 
here, therefore, are preliminary and are primarily qualitative in nature. Beyond 
the wing tip, those inner grid points which define the airfoil section, coalesce to 
essentially a single location as shown in Figure 5.30. This line which extends from 
the wing tip to the outer boundary in the spanwise direction is termed the “singular 
line”. Boundary conditions must be imposed at all the coalesced points located on 
the singular line. All points immediately surrounding the singular line are used to 
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Figure 5.29: View in the plane of the rotor of the 3-D, two-zone grid used for com 
putations of wing/rotor interaction. 




Figure 5.30: View of a vertical x — z plane of the grid outboard of the wing tip of 
the 3-D, two-zone grid used for computations of wing/rotor interaction. 


determine the average value of each of the flow properties. This is used to specify 
the boundary conditions on the singular line. This region, which is of no major 
concern in the freestream flow solutions, is seen to cause problems in the wing/rotor 
interaction calculations. This will be discussed in greater detail later. 

For these calculations, a uniform A p is again imposed across the rotor, but, for 
this case, a more representative value of O.Olpoo is used. 

Figure 5.31 shows the Mach number contours in a vertical x — z plane at the mid 
semi-span location. For this calculation, the Mach contours pass smoothly through 



Figure 5.31: Mach number contours in a vertical plane through the mid semi-span 
location for 3-D wing/rotor interaction, where A p = O.Olpoo and the rotor diameter 
is 4.78. 

the rotor. Figure 5.32 is the corresponding plot for the static pressure contours. 

Figures 5.33 to 5.35 are all views of flow properties in a vertical plane that runs 
spanwise through the wing mid-chord. The rotor location is superimposed on these 
plots. Intuitively, the flow above the wing seems to have been computed correctly. 
The flow stagnates in the mid-chord region at the semi-span location corresponding 
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Figure 5.32: Pressure contours in a vertical plane through the mid semi-span location 
for 3-D wing/rotor interaction, where A p = O.Olpoo and the rotor diameter is 4.78. 

to about y/c = 2.0. From this point the flow is divided into two opposite spanwise 
directions - towards the tip and towards the wing root. 

Beyond the wing tip, however, the singular line of the grid is affecting the results. 
It seems to be acting as a source, producing a discontinuity in velocity along its 
length. The reason for this behavior is as yet unknown. Figure 5.36 is a view 
of the Mach number contours in a vertical x — z plane beyond the wing tip. A 
large gradient in flow velocity is produced in a very localized region. The singular 
line boundary conditions are treated identically for the 3-D freestream computations, 
and, as previously mentioned, no anomaly was found for that case (see Section 5.2.2). 
It is possible that there may be some inconsistency in the application of the boundary 
conditions for the rotor, and that this may be affecting the computation on the inner 
boundary that contains the singular line. Various attempts at correcting the situation 
have so far been unsuccessful. One way to eliminate the problem is to develop a grid 
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Figure 5.33: Velocity vectors in a vertical, spanwise plane through the win* 

mid-chord for 3-D wing/rotor interaction, where A p = O.Olpoo and the rotor di 
ameter is 4.7S. 
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Figure 5.34: Velocity vectors around the wing tip in a vertical, spanwise plane 
through the wing mid-chord for 3-D wing/rotor interaction, where A p = O.Olpoo 
and the rotor diameter is 4.78. 

with no singular line. Although grid generation becomes considerably more complex, 
this approach is currently being pursued. 

5.5 Wall Jet in Two Dimensions 

Some preliminary results for tangential blowing near the leading edge have also been 
obtained. The furthest forward slot on the 0.18 scale model soon to be tested at 
NASA Ames’ OARF is located at 8.1% of the airfoil chord from the leading edge. A 
tangential jet at this location is modeled by imposing a jet velocity at the appropriate 
grid locations. This jet velocity is derived by assuming isentropic expansion from the 
plenum within the airfoil where the total pressure is specified. Assuming a plenum 
pressure 8% above freestream ambient pressure, the jet velocity at the plenum exit 
is about M = 0.3. Figure 5.37 is a plot of the velocity vectors for a 2-D wing/rotor 
interaction case with tangential leading edge blowing. A uniform pressure rise of 
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Figure 5.35: Mach number contours in a vertical, spanwise plane through the wmg 
mid-chord for 3-D wing/rotor interaction, where A p = O.Olpoo and the rotor diameter 
4.78. 


Figure 5.36: Mach number contours in a vertical x — z plane beyond the wing tip for 
3-D wing/rotor interaction, where A p = O.Olpoo and the rotor diameter is 4.78. 








Figure 5.37: Velocity vectors for a 2-D wing/rotor interaction case with tangential 
wall blowing, where the plenum pressure is 1.08poo, A p = 0.05poo across the rotor, 
and Mqo = 0.001. 

0.05poo across the rotor is specified for this hover calculation. It can be seen that the 
magnitude of the jet velocity far exceeds any of the rotor-induced flow velocities. 

Figure 5.38 is a close-up view of the jet velocity profiles in the downstream of the 
jet exit. This case assumes a freestream flow of Mqo = 0.2 at —90°, in addition to a 
rotor pressure rise of 0.03poo. The plenum pressure, as before, is 1.08poo- Of interest 
in this plot is the evolution of the velocity profiles as the outer flow is entrained by 
the jet through momentum transfer due to viscous effects within the boundary layer. 
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Figure 5.38: Close-up view of velocity vectors for a 2-D wing/rotor interaction case 
with tangential wall blowing, where the plenum pressure is l.OSpoo, A p = O.OSpoo 
across the rotor, and M 0 0 = 0.2. 
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Chapter 6 

Conclusions and 
Recommendations 

6.1 Conclusions 

The research undertaken to date has demonstrated the feasibility of analyzing the 
tilt rotor flowfield using the thin-layer Navier-Stokes equations and an implicit, finite 
difference technique to solve them. Computational fluid dynamics promises to yield 
the most accurate, theoretical prediction of download yet. Good results have been 
obtained for the V-22 airfoil and wing in a low subsonic Mach number flow ( M ^ = 
0.2) at -90 degrees. The use of multiple grid zones has enabled a more accurate 
representation of the full rotor flowfield. A rotor alone in two and three dimensions 
has been modeled using two zones. The results are very good and can be further 
improved by grid refinement. Wing/rotor interaction has been modeled in such a way 
that rotor and wing flows are partially coupled. Although the A p across the rotor 
remains fixed, the pressures in the rotor plane are influenced by the wing which lies 
less than a chord length below the rotor. Good results in two dimensions have been 
obtained for a lifting rotor in hover. In three dimensions, the flow computed over the 
wing looks intuitively correct, but beyond the wing tip, the presence of the singular 
line in the grid causes a local discontinuity in the computed flowfield. This problem 
is currently being investigated. Finally, preliminary results obtained for tangential 
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leading edge blowing in two dimensions look very encouraging. 

6.2 Recommendations 

6.2.1 Near Term 

It is first desired to improve the three-dimensional wing/ rotor interaction solution by 
solving the problem associated with the singular line in the grid outboard of the wing 
tip. Perhaps an inconsistency in the boundary conditions exists that could be alle- 
viated by improved implementation. Alternatively, a new, albeit more complicated, 
grid could be generated which has no singular line. 

The elliptic grid solver should be employed with the boundary control suggested 
by Sorrenson and Steger. This would smooth the interior grid point distributions 
and remove the kinks, thereby improving spatial accuracy of the solution. 

Modification of the code to allow for non-uniform, radial distributions of pressure 
rise and swirl is a fairly straight-forward task. 

A turbulence model should be implemented. A Baldwin-Lomax model would 
be adequate for the boundary layer, and Roberts’ model would be suitable for the 
curved jet region. 

Tangential leading edge blowing over the finite wing should also be implemented. 

The above modifications should be followed by detailed flow calculations and 
comparisons with available experimental results obtained in scale model tilt rotor 
tests by the Rotorcraft Aeromechanics Branch of the Full Scale Aerodynamics Divi- 
sion of NASA Ames Research Center. Effect of rotor thrust coefficient, swirl velocity, 
swirl direction, jet blowing strength and location, on the wing download should be 
studied. 

It is the intention of the authors to pursue these near term recommendations 
during the coming year. 
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6.2.2 Far Term 


In the far term, there are several improvements to the methodology that could be 
introduced to yield a computation that is more accurate and more representative of 
an actual tilt rotor vehicle. 

One of first areas of focus would be the rotor. Eventually, at some point in the 
future when computational resources become available, it would be desirable to per- 
form a time-accurate Navier-Stokes calculation about each of the rotating blades of 
the rotor and to couple that with the solution of the rest of the tilt rotor flowfield. 
In the interim, however, incorporating the momentum theory/blade element anal- 
ysis into the Navier-Stokes flow solver would be a considerable improvement. This 
would allow radial and azimuthal variations of flow properties on the rotor disk that 
would also be influenced by the wing. A p, for example, would not be fixed prior to 
the Navier-Stokes calculation but would result from it. The improved rotor model 
could be incorporated into the boundary conditions of the flow solver at a common 
boundary between two zones. Alternatively, by adding momentum source terms to 
the momentum equations of the Navier-Stokes equations, the rotor could be defined 
at any desired set of grid points in the mesh interior, thus eliminating the need for 
a zonal approach to rotor modeling. It is recommended that both of the approaches 
be attempted, and the results compared. 

Modeling the engine pylon and the fuselage are obvious improvements that would 
help yield a more accurate representation of the actual tilt rotor aircraft. This, of 
course, would entail complex grid generation. In the same vein, incorporation of a 
trailing edge flap in the airfoil grid would be helpful in the analysis of the effect of 
download due to flap deflection. 

A greater degree of accuracy in the modeling of the separated flow region would be 
attained if the full Navier-Stokes equations are solved. Also, an improved turbulence 
model developed specifically for bluff body flows would definitely improve the results 
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even further. A computation which is second-order accurate in time might also be 
attempted to better resolve the unsteadiness of the flow. 

Implementation of the above recommendations is a challenging task, but one that 
would be invaluable in the pursuit of the complete understanding of the complicated 
tilt rotor flowfield. 
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